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1.  Executive  summary  of  report 

This  report  summarizes  the  work  done  at  the  National  Bureau  of 
Standards  (NBS)  under  its  project  WHERE  in  support  of  the  US  Army  Advanced 
Materiel  Concepts  Agency's  (USAAMCA)  project  MNPLS  (Micro  Navigation  and 
Position  Location  System) . 

MNPLS  is  a system  that  utilizes  a time-division  multiplexed  frequency 
that  is  shared  by  n units , all  of  which  are  synchronized  by  a suitable 
electromagnetic  signal  during  each  time  interval  AT  . Each  unit  is  assigned 
one  or  more  time  slots  within  AT  during  which  it  emits  an  electromagnetic 
signal.  This  signal  in  turn  is  received  by  several  but  not  necessarily  all 
other  units  in  the  field  which  measure  the  time  of  flight  of  the  EM  signal; 
from  the  knowledge  of  the  time  and  the  assignment  of  the  time  slots  for  each 
unit,  the  distance  between  the  sending  unit  and  the  receiving  unit  can  be 
inferred . 

For  NBS  the  problem  to  be  solved  was  twofold: 

(1)  To  determine  the  feasibility  of  such  a system  in  the  presence  of 
measurement  errors.  If  all  measurements  were  exact  three  units  in  the  plane 
or  4 units  in  three-space  would  suffice  to  locate  another  (or  all  other)units. 
In  the  presence  of  error  more  measurements  are  needed  and  the  "best"  solution 
that  minimizes  total  error  in  some  mathematical  sense  has  to  be  found.  The 
fundamental  problem  for  NBS  was  to  show  that  the  overall  error  does  not 
increase  unduly  as  the  system  develops  in  time. 
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(2)  To  develop  a computer  simulation  for  such  a system  that  allows  the 
investigation  of  real-time  scenarios  with  a sufficiently  large  number  of 
units  undergoing  battle-field  manoeuvres,  each  with  prescribed  movements 
over  suitably  long  times.  The  simulation  program  determines  the  distance 
between  any  two  units  i and  j , perturbs  that  distance  according  to 
suitable  assumptions  about  the  measurement  errors,  and  supplies  to  the 
master  computer  (the  program  that  determines  positions  from  the  distance 
measurements)  the  distances  as  though  they  came  from  actual  field  measure- 
ments . 

The  output  of  the  program  consists  of  the  successive  locations  of  all 
units,  the  errors  between  true  position  and  calculated  position,  and  a 
measure  of  confidence  derived  from  the  available  measurements  but  not  from 
the  knowledge  of  exact  position  (which  is  only  available  to  the  simulation 
program  but  not  to  the  position  location  algorithm) „ 


NBS  has  designed  a set  of  mathematical  algorithms  for  position  location 
determination  that  allow  the  determination  of  successive  positions  of  n 
field  units,  from  inputs  containing  only  distances  between  units,' or  perhaps 
those  distances  accompanied  by  independent  measurements  of  differences  in 
height  above  sea  level*  The  experiments  based  on  actual  deployments  in  the 
Boston  area  show  that,  with  a suitable  number  of  fixed  units,  the  average 
error  of  position  overall  units  can  be  kept  well  within  acceptable  limits. 
Based  on  a range-error  standard  derivation  of  6m  , and  using  all  83  units 
as  reporters,  the  maximum  error  during  2 hours  of  actual  time  remained  under 
110m.  The  average  error  was  4.5m  over  all  units,  with  only  3 units  "lost" 
more  than  8 times  out  of  240  tries.  When  only  the  39  fixed  units  were  used 
as  reporters,  the  average  error  was  4.7m,  the  maximum  was  40m,  but  15  units 
were  "lost"  more  than  50  times  out  of  240.  (The  term  "locators"  refers  to 
those  users  that  do  actually  receive  a signal  from  the  locatee,  and  do 
report  a time -of -flight  to  the  master  computer;  the  term  "reporters"  refers 
to  those  users  that  have  been  instructed  to  report  time-of -flight  measure- 
ments to  the  computer.  Because  of  "line-of -sight"  requirements,  not  all  of 
the  latter  may  be  able  to  serve  as  locators  for  any  given  location-operation.) 

A simulation  package  WHERSM  was  designed  that  simulates  the  operation  of 
a Micro-Navigation  and  Position  Location  System.  This  program  allows  one  to 
stage  a wide  variety  of  movements  for  any  battle-field  scenario;  it  generates 
the  inter-unit  distances,  perturbs  the  ranges  according  to  agreed  upon  error 
distributions  and  provides  the  links  to  the  position  location  algorithms 
which  were  also  developed  at  NBS  for  this  purpose. 
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The  position  location  algorithms  fall  into  three  distinct  categories: 


(1)  Two-dimensional  algorithms.  Where  the  terrain  is  flat,  no  altitude 
reference  is  necessary  and  the  resulting  algorithms  become  two-dimensional. 

2 2 2 

The  sura  £w.(r,  - d.  ) is  minimized,  where  the  sum  extends  over  some  or 
i i l 

all  reporting  units  and  the  weights  are  determined  by  the  estimates  of 
variance  of  each  of  the  reported  distances  d^ 

(2)  Pseudo  three-dimensional  algorithms.  Where  the  altitudes  are  so  small 
in  relation  to  the  horizontal  distances  that  the  slant  distances  are  nearly 
the  same  as  the  horizontal  distances  it  becomes  necessary  to  reduce  the 
position  location  problem  to  the  two-dimensional  one  above.  Not  only  is  the 
three-dimensional  problem  ill  defined  mathematically,  but  we  have  also  a 
genuine  multiplicity  of  solutions  if  all  locators  are  nearly  in  a plane  and 
trying  to  locate  a unit  outside  the  plane.  The  determination  of  heights  that 
are  small  against  the  horizontal  distance  is  impossible  within  any  reasonable 
limits  of  accuracy  and  the  altitude  must  therefore  be  determined  either  by 
barometric  reference  or  a reference  to  a digitized  terrain  map.  For  example, 
at  a distance  of  1km,  a range  uncertainty  of  only  6m  leads  to  a height 
uncertainity  of  110m. 

(3)  Fully  three-dimensional  algorithms  are  being  used  when  the  altitudes  are 

comparable  to  the  horizontal  distances.  Again,  a penalty-function  approach 

2 2 2 

is  used  that  minimizes  £ w^(d^  " ri  ) with  suitably  determined  weights 
depending  on  the  variance  of  the  measurement  and  the  sum  being  taken  over 
some  or  all  reporting  units. 
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The  following  results  were  obtained: 


o The  test  runs  conducted  on  the  Boston  deployment  show  that  stability 
of  the  position-location  process  is  indeed  achievable  as  long  as  some  units 
remain  fixed. 

o Computer  size  for  position  location: 

Runs  were  all  carried  out  on  the  UNIVAC  1108  at  NBS . For  this  machine, 
the  formula  for  computing  the  number  of  computer  words  required  by  the 
position-location  algorithms  as  they  are  currently  programmed  is  as  follows: 

26N  + 7R  + 3874,  • 

where  N » number  of  units  in  the  field,  and  R = the  maximum  number  of 
ranges  being  reported. 

As  a matter  of  possible  interest,  the  number  of  additional  words 
required  by  the  simulation  program  (into  which  the  position-location 
algorithms  were  inserted)  is: 

9 IN  + Ng  + 4002, 

where  N = the  number  of  subcycles  in  a cycle  of  the  MNPLS  (see  Working 

u 

Paper  No.  9)  „ Furthermore,-  the  preprocessor  to  the  simulation  requires 
39,714  computer  words. 

o Timing  for  position  location  operation: 

Several  runs  are  available  from  which  the  time  required  by  the 
algorithms  can  be  estimated.  One  run  was  done  using  LSS(3d),  with  20  locators 
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it  required  approximately  16«5ms  per  location  operation.  A run  with  15 
locators  using  slant  range  reduction  took  8-12ms  ; with  20  locators,  the 
computation  took  14  to  18ms.  For  the  Boston  deployment,  using  39  reporters 
(but  probably  less  than  15  locators  on  the  average  because  of  intervisibility 
restrictions),  8 to  14ms  were  required;  using  83  reporters,  17  to  25ms  were 
required.  Since  the  algorithms  are  not  optimized  for  running  time,  it  is 
reasonable  to  expect  some  decrease  in  these  times  for  the  operational  system. 

2 . Description  of  general  approach  used  in  this  report, 
a)  Simulation. 

The  simulation  of  the  battlefield  scenario  is  done  in  our  NBS-developed 
computer  program  WHERSM.  This  program  provides  for  a "real  life"  framework 
for  MNPLS  through  its  capability  for  moving  all  field  units  along  prescribed 
paths,  or  subjecting  these  to  random  changes  in  their  direction  of  movement. 
WHERSM  generates  all  the  inter-unit  distances (ranges) ; checks,  if  necessary, 
for  intervisibility  (by  look  up  in  an  intervisibility  matrix)  ; disturbs  the 
exact  ranges  according  to  known  error  distribution  laws;  and  transmitts  these 
"measurements"  to  the  position  location  algorithms  that  we  developed  here  at 
NBS  and  which  are  described  below. 

The  simulation  package  also  provides  a facility  for  monitoring  the 
operation  of  such  a system  under  a variety  of  movement  scenarios;  and 
ultimately  it  provides  for  a variety  of  outputs  that  are  desired  for  both 
checking  out  the  feasibility  of  the  entire  system  and  what  might  be  called 
the  ultimate  output  that  would  be  useful  to  the  field  commander  who  is 
interested  in  the  whereabouts  of  his  units.-  For  details  of  WHERSM,  the 
reader  is  referred  to  working  paper  No.  9 in  the  appendix  of  this  report. 
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(2)  Algorithms . 


The  basic  LSS  algorithm  that  we  finally  used  after  numerous  experi- 
ments designed  to  eliminate  competing  ones  is  the  so-called  "least-squares 
squared"  algorithm,  which  minimizes  the  penalty  function 

^ ,.2  2 N2 
E = £ w.  (d.  - r . ) 

where  i is  taken  over  all  of  the  locators  and  w^  is  a weight  that 
increases  with  decreasing  variance  of  the  observation  , the  reported 
distance  of  the  locatee  from  the  locator  with  index  i . The  quantity  d^ 
is  the  calculated  distance  between  (x^,  y^,  z ) and  (x,y,z)  the  estimated 
position  of  the  locatee.  The  solution  (x,y,z)  is  obtained  iteratively 
from  the  equations  BE/dx  = 0 , dE/dy  = 0 , BE/dz  = 0 once  a starting 
point  is  known. 

As  long  as  a trajectory  is  known  for  a given  locatee,  the  initial  guess 
is  simply  obtained  from  extrapolation.  If  the  unit  is  lost,  however,  a set 
of  linear  equations  is  solved  (LSL  method).  Details  of  the  LSL  method  can 
be  found  in  NBS  report  10663:  First  Interim  Progress  Report  on  NBS  Project 

WHERE  in  support  of  USA  AMCA  Project  MNPLS. 

The  appendix  contains  a list  of  variables  used  in  the  simulation  code. 
Since  this  is  the  part  of  the  NBS  system  that  interfaces  with  the  outside 
world,  we  believe  that  those  are  the  only  variables  that  have  to  be  listed. 
The  variables  in  the  algorithm  parts  are  internal  and  are  thus  not  referenc- 
able  from  outside. 
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Review  of  analyses  to  date. 


3.1,  General . 

✓ As  originally  conveived,  the  testing  was  to  comprise  two  phases: 

initial  testing  to  choose  the  best  position-location  algorithm,  and  final 
verification  testing  to  determine  how  well  the  chosen  algorithm  would  work 
on  the  Boston  deployment.  This  neat  plan  did  not  fit  the  real-life  flow  of 
events,  however:  every  stone  turned  over  revealed  two  more  yet  unturned. 

The  process  of  developing  a new  system  is  bound  to  contain  such  loops.  For 
example,  one  does  not  know  how  to  optimize  the  algorithm  until  one  knows 
the  system  design,  but  the  system  design  is  sensitive  to  the  effectiveness 
of  the  chosen  algorithm.  Thus  the  two  phases  gradually  blurred  into  each 
other.  Furthermore,  in  the  limited  time  allotted  for  this  study,  it  seemed 
that  effort  would  be  better  spent  getting  reasonable  answers  to  more  of  the 
real  questions,  rather  than  "optimizing”  on  a smaller  subset  of  the  ques- 
tions while  assuming  particular  (and  questionable)  answers  to  the  others. 

For  these  reasons,  a strictly  chronological  account  of  our  analyses 
would  not  be  very  helpful.  Instead,  it  seems  reasonable  to  describe  the 
test  results  by  first  building  a framework  of  the  questions  to  be  answered, 
and  then  describing  the  tests  performed  and  the  answers  obtained.  The  next 
section  consists  of  a set  of  questions;  the  following  section  presents  the 
answers  to  date,  as  gleaned  from  the  experimental  program  and  from  certain 
theoretical  studies  undertaken  specifically  to  support  that  program. 

3.2,  Questions  pertinent  to  the  MNPL  system. 

The  general  question  that • constitutes  the  focus  of  a feasibility  study 
is:  Will  the  proposed  system  perform  satisfactorily  under  the  conditions  of 
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its  intended  use?  This  question  is  however  a very  broad  one;  it  asks 
essentially  whether  the  final  outcome  of  the  feasibility  study  is  "yes", 

"no"  or  "maybe” , and  thus  does  not  give  much  guidance  in  conducting  the 
study.  Thus  one  is  led  to  pose  narrower  queries: 

a.  General  Questions 

(1)  Under  highly  optimistic  conditions,  will  the  system  work?  (If  not, 
terminate  the  study.) 

(2)  Under  what  reasonable  near-minimum  conditions  can  the  system  be  made 
to  work? 

(3)  How  well  does  it  work? 

b . More  Specific  Questions 

The  next  question  leads  to  a whole  set  of  sub-questions,  which  are 
listed  along  with  the  main  one. 

(4)  What  position-location  algorithm(s)  should  be  used?  For  each  candidate 
algorithm,  consider: 

(a)  Computer  time  required? 

(b)  One-step  accuracy? 

(c)  Error  propagation  (sensitivity  to  locators1  position  errors)? 

(d)  Sensitivity  to  occasional  very  large  errors  in  the  range 
measurements? 

(e)  Affected  by  dimensionality? 

(f)  Need  an  initial  estimate  of  position? 

The  next  set  of  questions  refers  to  restrictions  placed  on  the  units  in 
the  system. 


(5)  How  many  (if  any)  units  must  remain  fixed? 

(6)  Must  they  be  permanently  fixed  in  known  positions? 

(7)  How  many  reporters  are  required,  to  keep  the  errors  small? 

(8)  Will  ground  units  have  to  carry  altimeters?  Instead,  will  a digitized 
map,  stored  in  the  computer,  suffice? 

(9)  How  accurate  must  altimeters  be  (both  ground  and  airborne)? 

Questions  (10)  - (19)  refer  to  design  of  the  algorithm. 

(10) How  often  should  units  report? 

(11) How  many  aircraft  are  required,  in  order  to  use  a three-dimensional 
algorithm  on  ground  units? 

(12) How  should  one  handle  near-planarity  (when  any  aircraft  present  are 
flying  very  low)? 

(13) When  should  the  algorithm? 

(a)  ignore  height  differences? 

(b)  reduce  slant  ranges  to  horizontal  distances? 

(c)  use  a 3-dimensional  method? 

(d)  use  3d  for  aircraft,  slant  range  reduction  for  ground  units? 

(14) How  should  the  locators  be  chosen? 

(15) How  weight  the  locators? 

(a)  What  is  the  distribution  of  range  errors? 

(b)  How  handle  the  occasional  occurrence  of  range  measurements  with 
extremely  large  errors? 

(c)  How  estimate  the  accuracy  of  an  estimated  position? 

(16) Should  the  locators  be  "relocated1'  ? (i.e.,  should  the  locators' 
positions  be  adjusted,  based  on  discrepancies  between  measured 
ranges  and  calculated  ranges  to  the  best  estimate  of-locatee's 
position?) 
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(17)  Use ' only  fixed  units  as  reporters?  How  soon  after  they  become  fixed? 

(18)  Can  the  algorithm  determine,  (satisfactorily)  which  units  are  fixed? 

If  not,  could  a moving/stopped  indicator  bit  be  included  in  the 
transmission? 

(19)  After  the  system  has  been  running  for  a while,  and  errors  have 

accumulated,  is  it  possible  to  ,,restartM  the  system i.e.  to  treat 

the  set  of  (currently)  fixed  units  in  toto,  simultaneously,  allowing 
all  position  coordinates  to  vary,  in  order  to  converge  to  good  estimates 
for  their  positions? 

3.3.  Conduct , progress  , and  results  of  the  test  program. 

The  first  few  questions,  above,  suggest  a beginning  emphasis  on 
favorable  conditions.  The  rationale  is  that,  when  very  little  is  known 
about  a system,  more  is  learned  faster  by  starting  with  conditions  where  the 
system  works  tolerably  well,  and  then  investigating  the  effects  of  departures 
from  these  conditions.  Accordingly,  the  first  major  project  was  to  devise, 
develop,  program,  and  test  several  candidate  algorithms  for  determining 
the  position  of  a locatee,  using  the  (perturbed)  ranges  to  a set  of  locators, 
and  assuming  the  locators’  positions  to  be  known  exactly.  That  is,  attention 
was  focussed  first  on  questions (4)  . Some  of  these  algorithms  were  suitable 
for  2 or  3 dimensions,  others  were  usable  in  2 dimensions  only.  Methods 
considered  were:  Linear  Method  (LM) , Smallest  Tangent  Circle  (STC),  Least 

Squares  Linear  (LSL) , Minmax,  Least  Squares(LS),  and  Least  Squares  Squared 
(LSS).  The  first  two  methods  use  exactly  3 locators  in  2 dimensions,  or  4 
locators  in  3 dimensions.  The  next  two  are  generalizations  of  the  first  two. 


// 


W I 


' 


respectively , to  an  arbitrary  number  of  locators.  Least  squares  and  LSS 
are  penalty-function  methods.  All  are  fully  described  in  earlier  sections, 
in  Working  Papers  Nos.  7 and  8,  and/or  in  NBS  Report  10663,  First  Interim 
Progress  Report  on  Project  WHERE  in  support  of  USAAMCA  Project  MNPLS . 

Tests  were  first  conducted  in  two  dimensions.  STC  turned  out  consider- 
ably more  accurate  than  LM,  but  slower;  LS  and  LSS  agreed  closely,  and 
were  more  accurate  than  any  others,  but  slower  than  all  but  Minmax9  Since 
Minmax  was  no  more  accurate,  and  several  times  more  time-consuming,  than 
LSS,  it  was  dropped  from  further  consideration.  Also,  since  LM  was  inferior 
to  STC,  it  was  dropped.  Thus  the  survivors  at  this  point  were:  STC,  LSL,  LS, 
and  LSS. 

The  next  sequence  of  tests  was  directed  to  question  (4)(c)  . For  these 
runs,  errors  were  applied  to  the  positions  of  the  locators,  corresponding 
to  a point  in  time-  after  the  system  has  been  operating  long  enough  to 
accumulate  such  errors.  In  these  runs,  STC  fell  down  markedly.  LS  and  LSS 
did  better  than  LSL,  but  the  latter  was  kept  since  it  alone  (of  the  three 
methods  still  in  contention)  does  not  need  a starting  value  — i.e.,  is  not 
an  iterative  method  — and  therefore  can  provide  a position  estimate  when 
the  locatee  first  joins  the  net  or  has  been  out  of  touch  for  some  time. 

Three-dimensional  tests  confirmed  the  pattern  of  results  above,  and 
also  accented  the  importance  of  geometry:  when  the  locators  are  nearly 
coplanar  (collinear  in  2 dimensions),  LSL  often  produces  a very  poor  estimate, 
because  the  planes  (lines)  that  it  tries  to  fit  are  nearly  parallel.  LS  and 
LSS  also  have  problems,  but  not  as  severe  as  LSL.  These  problems  were 
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accentuated  in  3d,  simply  because  locators  on  the  surface  of  the  earth  are 
relatively  coplanar;  the  corresponding  situation  in  2d  would  have  the  locators 
strung  out  nearly  on  a single  line,  a less  likely  situation. 

In  terms  of  computer  time  LSL  was  best,  with  LSS  next.  Since  LSS  was 
devised  to  perform  like  LS  but  faster,  these  two  were  then  compared  for 
accuracy.  72  trials  were  run  on  each  method,  with  18  users  scattered  over 
.an  8km  by  16km  area,  with  altitudes  to  9km,  for  each  of  two  different  values 
of  a (the  standard  deviation  of  the  range  error).  Each  of  the  18  units 
served  as  locatee  4 times  for  each  method  and  each  value  of  a . The  results 
were  as  follows;  For  a = lm  , the  position  estimates  were  in  error  by  as 
much  as  3.32m  , but  the  difference  between  the  LS  estimate  and  the  LSS 
estimate  was  never  more  than  0.0004m  „ For  a * 10m  , the  errors  reached 
31. 9m,  but  the  two  methods  gave  estimates  that  differed  by  no  more  than  0.045m. 
As  a result  of  this  test,  it  was  concluded  that  LS  and  LSS  could  be  used  inter- 
changeably so  far  as  accuracy  is  concerned.  Since  LSS  was  considerably 
faster,  LS  was  dropped.  The  final  result  thus  was;  Use  LSS;  if  no  starting 
value  (initial  estimate  of  position)  is  available,  use  LSL  to  obtain  one. 
(Complications  develop  when  there  is  a choice  between  2-  and  3-dimensional 
versions,  as  we  will  see.) 

The  only  part  of  Question  (4)  yet  unanswered  is  (4)(d):  sensitivity  to 
occasional  large  errors  ( outliers)  in  the  measured  ranges.  This  is  a 
severe  problem  for  STC  and  Minmax,  since  STC  uses  only  three  ranges  and 
Minmax  concentrates  on  minimizing  the  maximum  discrepancy.  For  LSL,  LS,  and 
LSS,  the  matter  is  treated  in  a statistical  sense,  as  is  commonly  done  in 
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fitting  regression  curves.  The  true  position  is  estimated,  and  then  the 
discrepancies  between  measured  and  calculated  ranges  are  examined.  Any 
discrepancies  significantly  larger  than  the  average  are  attributed  to  faulty 
measurements.  These  measurements  are  ignored,  and  the  cycle  repeated. 

This  approach  is  sure  to  detect  all  really  bad  errors,  provided  they  occur 
infrequently,  unless  there  is  exactly  the  minimum  number  of  ranges  necessary 
to  produce  an  estimate.  It  was  been  shown  to  be  quite  efficient,  in  the 
sense  that  the  occasional  rejection  of  measurements  "apparently"  but  not  really 
contaminated  with  large  extraneous  errors  does  not  appreciably  degrade  the 
accuracy  of  the  estimation  procedure.  The  actual  cutoff  point  for  rejection 
is  a variable  that  will  need  further  investigation  as  part  of  a system  develop- 
ment effort. 

Questions  5,6,  and  7 pertain  to  the  set  of  reporters.  (The  term 
"locators"  refers  to  those  users  that  do  actually  recieve  a signal  from  the 
locatee,  and  do  report  a time-of-f light  to  the  master  computer;  the  term 
"reporters"  refers  to  those  users  that  have  been  instructed  to  report  time- 
of-f  light  measurements  to  the  computer.  Due  to  the  requirement  for  a "line- 
of-sight",  not  all  of  the  latter  may  be  able  to  serve  as  locators  for  any 
given  location-operation.)  In  order  to  obtain  broad(rather  than  situation- 
specific)  answers  to  these  questions,  it  was  necessary  to  set  up  a number  of 
different  scenarios,  with  different  numbers  of  users  and  different  propor- 
tions stopped.  This  diversity  was  created  by  employing  several  kinds  of 
"randomness"  . Users  were  assigned  to  randomly  chosen  starting  points  within 
a specified  rectangle,  with  aircraft  heights  specified  or  randomly  chosen; 
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sequences  of  azimuths  were  chosen  randomly  over  a 90°  range , and  azimuth 
change  times  were  drawn  from  an  exponential  distribution.  Stopping  and 
starting  times  were  drawn  from  exponential  distributions,  chosen  to  achieve 
specified  values  for  the  expected  (long-term  average)  proportion  of  users 
stopped  at  any  given  instant,  and  for  the  mean  times  spent  in  ’'fixed" 
status  and  in  "moving" . Intervisibilities  were  similarly  randomized,  to 
achieve  specified  values  for  the  average  time  two  users  maintain  communication 
and  the  average  time  they  remain  out  of  Communication.  These  features  were 
used  in  various  combinations  to  investigate  most  of  the  remaining  questions. 

The  results  of  those  investigations  are  presented  below. 

If  all  the  reporters  are  moving,  strange  things  can  happen.  Small 
(random)  errors  in  locating  a given  reporter  propagate  into  the  location  of 
other  reporters,  so  that  each  remains  well-located  relative  to  the  others 
but  the  whole  cluster  is  displaced(and  perhaps  rotated)  far  from  its  true 
position.  This  effect  is  reasonably  well  understood,  and  has  been  clearly 
demonstrated.  One  case  involved  30  users,  moving  more  or  less  in  the  same 
direction.  12  users  were  reporting,  with  perfect  intervisibility  (2  dimensions) 
By  the  time  100  location  operations  had  been  done  on  each,  they  had  moved 
0.9  to  4.1km  . The  estimated  position  of  the  cluster  as  a whole  was  3.3km 
from  its  true  location,  but  each  element  was  located  within  150m  of  its  true 
position  relative  to  the  others.  (Range  errors  were  less  than  2m  in  this 
case.)  The  problem  can  be  seen  clearly  in  the  following  simple  example:  If 
all  of  a collection  of  fixed  users  suddenly  begin  to  move,  each  with  the 
same  velocity  vector,  or  all  of  a collection  of  users  moving  with  a common 
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velocity  vector  suddenly  stop,  no  algorithm  based  on  inter-unit  ranges  can 
see  the  change,  because  the  inter-unit  ranges  don't  change.  Therefore,  the 
algorithm  will  leave  the  users  fixed  in  the  first  case,  and  will  continue  to 
move  them  in  the  second.  This  is  a special  case,  of  course,  but  the  same 
principle  applies  to  the  average  velocity  vector  when  the  units  begin  to 
move  with  different  velocities  - or  more  generally,  jto  the  average  change  in 
the  velocity  vector  when  they  (more  or  less)  all  change  direction  at  about 
the  same  time. 

These  results  made  it  plain  that  one  needs  some  users  fixed.  The  next 
question  was,  how  many? 

Several  runs  were  made  with  different  sets  of  fixed  users,  under  two 
conditions:  first,  that  the  computer  does  not  know  which  users  are  fixed, 
and  second,  that  it  does  know(and  can  use)  this  information.  The  first  run, 
with  range  errors  distributed  uniformly  from  -3m  to  +3m  , involved  15 
motorized  and  15  foot  units;  in  each  set  of  15,  5 were  fixed  and  reporting,  5 
were  moving  and  reporting,  and  the  remaining  5 were  moving  but  not  reporting. 
After  220  cycles  (109  minutes),  the  xy  errors  ranged  from  500m  to  10,240m, 
and  were  distributed  into  all  four  quadrants  (more  specifically,  14  users 
were  placed  NW  of  their  true  locations,  5 SW,10  SE,  and  1 NE  ) . For  the 
second  run,  the  subcycle  At  was  reduced  by  a factor  of  5,  to  see  if  locating 
the  users  more  often  would  improve  the  accuracy.  After  613  cycles  (60  min), 
the  errors  ranged  from  570m  to  2200Qnu  (The  errors  grew  more  slowly  per 
cycle,  but  faster  in  real  time.)  Next,  4 of  the  fixed  users  were  assumed 
"known"  - i.e.,  no  estimation  of  their  positions  was  allowed  - to  provide  an 


anchor  o£  sorts.  After  239  cycles  (116  min.),  the  average  error  over  all 
users  was  9.4m  , and  the  maximum,  was  23m  ; over  all  cycles,  the  average  was 
7.2m  , the  maximum  82m  * But  there  is  some  evidence,  that  the  errors  fluctuate 
considerably  with  time:  in  the  set  of  ten  cycles  prior  to  the  above  (i.e., 
cycles  221-230) , the  average  was  23m  ; the  maximum  69m  ! 

At  this  point,  relocation  (Question  16)was  tried.  (Refer  to  Working 
Paper  No.  8 for  technical  details.)  Tentative  conclusions  from  the  dozen  or 
so  runs  are:  In  the  beginning,  relocation  increases  the  errors.  Eventually, 

the  errors  without  relocation  become  larger  than  those  with  relocation:  how- 
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ever,  this  did  not  show  up  until  about  2 hours  of  real  time  had  been 
simulated,  by  which  time  the  fast-moving  users  had  separated  from  the  slow- 
moving  and  fixed  users  by  about  20km  . (At  this  time,  the  overall  average 
and  maximum  errors  were  19m  and  557m  with  relocation,  and  34m  and  1036m  with- 
out. Least  Squares  was  also  tried;  the  average  errors  were  comparable,  but 
the  maxima  were  570  and  1695  respectively.) 

Perhaps  more  instructive  were  the  shorter  runs  with  differing  values  of 
the  fraction  mentioned  in  Working  Paper  No.  8,  which  measures  the  extent  to 
which  relocation  is  carried.  Runs  were  done  with  this  fraction  set  at  0„5, 
0.2,  0.1,  and  0 (i.e.,  no  relocation).  The  results  showed  steady  improve- 
ment as  the  fraction  went  down,  with  smallest  errors  occurring  when  no  reloca- 
tion was  done.  This  shows  clearly  that  there  is  no  advantage  to  relocation 
for  the  first  100  cycles  (50  minutes)  at  least.  Still  open  (and  reasonable) 
is  the  possibility  that  relocation  should  be  instituted  for  a time,  now  and 
then,  with  a small  fraction,  to  help  remove  accumulated  errors. 


n 


However,  if  there  are  sufficiently  many  fixed  reporters,  one  can  probably 
do  better  by  periodically  going  into  a survey  mode  wherein  all  the  fixed 
units’  positions  are  simultaneously  adjusted,  effectively  ’’restarting”  the 
system. 

At  this  point,  it  seemed  clear  that  relocation  was  not  worth  its  price 
(more  than  doubling  the  algorithm’s  computer  time).  This  conclusion  was 
checked  once  more,  in  the  next  set  of  runs,  under  more  realistic  conditions, 
and  was  supported  by  the  results:  the  errors  were  smaller  (at  100  cycles) 
without  relocation. 

Parallel  efforts  by  ECOM  had  meanwhile  produced  a ’’best  engineering 

judgment"  range-error  distribution,  which  can  be  adequately  modeled  as 
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"Gaussian  with  variance  36  m ; mean  0 , except  that  a random  1%  of  measure 
ments  are  inflated  by  about  7 meters."  There  will  also  be  occasional 
extremely  large  errors;  a simple  preprocessor  can  weed  most  of  these  out  in 
the  field  system,  by  simply  checking  whether  the  range  measurement  is 
compatible  with  the  last  known  position  and  the  velocity  potential  of  the 
user.  This  error  distribution  was  implemented  in  the  balance  of  our  tests. 
Also  implemented  were  moving  boundaries,  which  can  be  set  to  advance  at  a 
speed  appropriate  for  foot-soldiers;  any  motorized  unit  which  reaches  the 
boundary  will  be  reflected  off  it,  so  that  the  motorized  units  do  not  stray 
far  from  the  foot  units.  Additional  features  implemented  at  this  time  were: 
a moving/stopped  indicator;  calculation  of  weights  which  reflect  the 
estimated  accuracy  of  a user’s  estimated  position;  the  ability  to  choose  as 
locators  those  users  with  greatest  weights,  or  those  which  have  been  stopped 
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for  a certain  number  of  location-operations;  improvement  of  estimates  for  the 
positions  of  fixed  users , by  combining  successive  estimates;  a functional 
representation  of  terrain;  use  of  an  independent  measure  of  the  z-coordinate 
(the  height),  with  random  errors;  algorithms  to  reduce  slant-range  distances 
to  planar  (xy-plane)  distances,  using  the  independent  measure  of  height 
difference  between  locator  and  locatee;  and  tests  to  determine  whether  a 
particular  location-operation  should  be  treated  by  slant-range  reduction  or 
by  a 3 -dimensional  method.  Again,  these  features  were  used  in  various 
combinations  for  the  runs  recounted  below. 

Before  these  remaining  runs  are  described,  the  rational  used  to 
"weight"  the  various  range  measurements  will  be  explained.  No  claim  is  made 
that  these  weights  are  optimal,  but  it  is  felt  that  they  are  reasonably  close 
to  the  best  values,  which  is  sufficient  for  a feasibility  study.  Suppose  half 
the  locators  had  the  same  x-coordinate  as  the  locatee,  and  half  had  the  same 
y-coordinate . Then  (assuming  ranges  to  be  large  relative  to  range  errors) the 
first  half  of  the  locators  is  useless  in  determining  the  x-coordinate,  and  the 
second  half  is  useless  in  determining  the  y-coordinate.  Thus  if  there  are  n 
locators,  n/2  are  used  for  each  coordinate,  and  so  we  have  n/2  estimates 
(assumed  independent)  of  each  of  the  two  coordinates  of  the  locatee* s 
position.  Consider  for  definiteness  the  x-coordinate.  If  the  locators'  posi- 
tion errors  are  unbiased  (i.e.,  have  mean  value  equal  to  zero),  the  best 
estimate  of  the  locatee* s x-coordinate  is  a weighted  average  of  these  n/2 
values:  weighted  by  the  reciprocals  of  the  variances  of  these  observations. 
Since  the  observational  error  is  the  sum  of  the  locator's  x-coordinate  error 
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and  the  range  measurement  error , its  variance  is  the  sum  of  the  position 

error  variance  (determined  when  the  locator  was  last  located)  and  the 
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measurement  error  variance  ( **  36  m ) 0 Thus  one  can  calculate  an  estimate 
of  the  x-coordinate  error  variance,  which  in  turn  is  used  to  calculate  the 
weight  for  this  user  when  he  serves  as  a locator. 

[ There  are  two  considerations  which  should  be  mentioned  here.  First, 
locators  are  not  in  general  split  so  nicely  along  perpendicular  lines.  Does 
this  matter?  Of  course  it  does,  in  the  following  sense:  if  more  locators 
serve  to  estimate  x than  y,  then  x is  likely  to  be  better  known  - i.e., 
to  have  smaller  variance  — and  y will  be  less  well  known.  On  the  average, 
however,  things  even  out:  if  azimuths  (from  locatee  to  locators) are  uniformly 

distributed  from  0 to  360*  , the  average  (or  expected)  variance  in  a given 
direction  turns  out  to  be  the  value  derived  above.  The  second  consideration 
is:  What  if  the  estimated  variances  for  the  locators  are  wrong?  Two 

consequences  follow:  (a)  To  the  extent  that  the  variances  are  different 
multiples  of  the  assumed  values,  the  estimate  of  position  is  less  than 
optimum  (because  the  true  relative  variances  should  be  used  to  get  the 
estimate) $ and  to  the  extent  that  the  variances  are  (as  a group)  larger  (or 
smaller)  than  assumed,  the  estimated  variance  of  the  locatee* s position  will 
be  too  small  (or  large,  respectively).  Point  (a)  is  not  likely  to  be  important 
since  it  would  take  very  large  discrepancies  to  affect  the  estimate  noticeably. 
Point  (b),  on  the  other  hand,  should  at  least  be  investigated.  One  technique 
is  to  act  as  though  the  true  variances  are  proportional  to  the  estimated  values 
with  an  unknown  constant  of  proportionality  (say  c)  . Then  c can  be 
estimated,  and  used  to  produce  a fair  variance  estimate.  Specifically,  if  each 
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range  observation  has  variance  co-^  (a^  *=  location  variance  + range  measure- 
ment variance),  then  the  best  estimate  for  the  locatee  can  be  found  without 
considering  c.  If  d^  represents  the  discrepancy  between  the  calculated 
range  from  this  position  to  the  assumed  position  of  the  i-th  locator  and 
the  corresponding  measured  range,  then  [ £ d^/a^  ^ t*ie  (multiplica- 

tive) correction  to  be  applied  to  the  variance  estimate  already  given.  (If 
c - 1 , this  value  should  also  be  ~1,)  This  is  a recent  result  and  the 
corresponding  modification  has  not  yet  been  made  to  the  algorithms.] 

A new  series  of  2-dimensional  runs  was  done,  with  30  users.  An  initial 
set  of  fixed  users  (reporters)  was  specified.  From  then  on,  units  became 
reporters  as  soon  as  they  had  been  fixed  for  ten  location— operations , and 
ceased  to  be  reporters  as  soon  as  they  began  to  move.  Average  times  moving 
and  fixed  were  set  at  15  min.,  so  that  about  15  (half  the  units)  were  stopped 
at  any  one  time.  A 2%-hour  (300  cycle)  run  was  done,  and  resulted  in  an  over- 
all average  error  of  3.2m,  and  a maximum  error  of  35raj  in  the  last  set  of  5 
cycles,  the  average  was  3.5m,  the  maximum  11m  . Next,  a run  using  longer 
cycles  (90  sec.  instead  of  30)  produced  errors  of  3.8  and  50.5  overall,  5.7 
and  33.8  for  the  last  5 cycles.  Next,  the  average  times  fixed  and  moving 
were  set  at  45  min.  and  this  run  repeated  (because  with  15  min.  times  and 
\\  min.  cycles,  many  units  didn’t  stay  stopped  long  enough  to  be  much  use  as 
locators).  The  error  figures  then  came  back  down  toward  those  of  the  first 
run:  3.4  and  33.3  overall,  5.9  and  22.5  for  the  last  5 cycles. 

At  this  point,  the  cycle  time  was  returned  to  30  sec.,  aid  a relocation 
run  was  tried,  with  fraction  0.2.  It  did  poorly:  errors  5.8  average  and 


35  maximum  overall,  with  8.4  and  17.5  for  the  last  5 cycles.  All  errors 
were  in  one  quadrant,  indicating  that  some  instability  had  been  introduced. 

Next,  a run  with  fraction  .001  was  done.  Since  this  fraction  is  so  small, 
the  magnitudes  of  the  actual  relocations  will  be  essentially  negligible; 
thus  the  effect  of  relocation  in  this  run  is  simply  to  allow  more  iterations 
for  the  algorithms.  (Remember  that  the  number  of  iterations  has  been  kept 
small  to  keep  the  running  time  down,  with  the  thought  that  range  errors  of 
the  order  of  6m  don't  justify  too  much  precision  in  locating  the  true  solution 
of  the  least  squares  problem.)  Only  100  cycles  were  done,  with  results  equal 
(up  to  roundoff  effects)  to  those  without  relocation.  This  supports  the 
conclusion  that  generally,  it  is  not  worth  doing  too  many  iterations.  (Of 
course,  how  many  is  "too  many"  depends  very  much  on  how  good  the  algorithm 
is*.) 

The  algorithm  was  modified  to  treat  aircraft  with  a 3-dimensional  method, 
while  continuing  to  treat  ground  units  2 -dimensionally . (This  is  possible 
since  only  ground  units  were  being  used  as  locators.)  Ten  aircraft  were 
added,  at  (random)  altitudes  ranging  from  128  to  2220m.  That  low-flying  aircraft 
presented  a problem  became  very  obvious,  since  the  algorithm  located  them  where 
they  belonged  at  some  times,  near  the  ground  at  others,  and  below  ground 
(mirror  image)  at  still  other  times.  One  more  attempt  was  made  in  this  series, 
using  only  ranges  that  were  at  most  10  times  height,  to  ensure  that  not  all 
lines  of  sight  were  flat.  This  helped  for  the  second-lowest  aircraft,  at  339m: 
it  had  an  average  (x,y)  error  of  19m,  with  a maximum  of  113m.  However,  the 
lowest  aircraft- generally  did  not  have  enough  acceptable (short  enough)  ranges 
to  estimate  its  location.  It  was  decided,  therefore,  to  do  a separate  systematic 


study  of  aircraft  location  errors.  This  will  be  discussed  below,.  The 
conclusion  is  already  clear , however , that  some  independent  height  informa- 
tion will  be  required  to  locate  low-flying  aircraft,  unless  other  aircraft 
are  serving  as  locators;  however,  the  use  of  aircraft  as  locators  will 
probably  require  very  stable  aircraft  and  good  tracking  techniques. 

A functional  representation  of  terrain  - i.e.,  z = f(x,y),  with  f(x,y) 
chosen  to  give  a reasonable  surface-was  implemented  at  this  stage,  to  test 
the  use  of  independent  pressure  references  (p.r.'s)  as  a source  of  height- 
difference  information.  After  consultation  with  AMCA  personnel,  these  were 
modeled  as  giving  height  with  an  error  that  had  zero  mean  and  a - 14m  (so 
that  the  difference  of  2 heights  will  have  c = 20m  ).  40  units  were  used, 
on  a terrain  that  varied  200m  in  elevation,  with  15  units  stopped  initially) 
again  only  units  that  had  been  stopped  for  10  location-operations  were  used 
as  locators.  The  x,y  coordinates  were  obtained  by  reducing  slant  ranges  to 
horizontal  distances  and  applying  2d  algorithms.  In  2%  hours,  the  overall 
average  error  was  3.5m,  the  maximum  29.2  ) for  the  last  5 cycles,  the 
average  was  4.8,  the  maximum  13.7  . Next  the  p.r.  a was  set  to  0 , with 
results  identical  to  another  run  where  the  terrain  was  flat  and  no  heights 
were  considered.  For  extraneous  reasons,  it  only  went  260  cycles  instead  of 
300)  it  had  errors  of  3.2  (avg),  28.2  (max)  overall)  3.3,  14.6  for  the  last 
5 cycles.  Next,  the  height  estimates  were  set  at  0 ; i.e.,  the  hills  and 
valleys  were  ignored  by  the  algorithm.  The  errors  were  slightly  smaller  than 
the  ordinary  run  using  p.r.'s:  3.35  and  28.2  overall,  4.4  and  14.3  for  the 
last  5 cycles.  This  is  not  surprising,  since  the  errors  introduced  by  the 
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p.r.'s  are  comparable  to  the  terrain  differences.  If  p.r.'s  don't  perform 
any  better,  perhaps  they  should  be  ignored  when  only  ground  units  are 
concerned.  However,  in  particular  situations,  the  terrain  heights  may 
introduce  systematic  errors  if  ignored,  which  may  then  propagate j the  p.r. 
errors  at  least,  have  the  advantage  of  being  random.  More  investigation  of 
this  possibility  should  be  undertaken  when  more  is  known  about  p.r.'s  . 

For  the  next  run,  terrain  roughness  was  tripled:  with  the  p.r.'s  the  errors 
were  essentially  the  same  as  on  the  flatter  terrain. 

Only  two  sets  of  runs  are  yet  to  be  described:  the  aircraft  error  study 
and  the  runs  on  the  Boston  terrain.  The  Boston  scenario  had  no  aircraft,  so 
it  was  treated  using  slant-range  reduction  with  p.r.'s.  Two  runs  using  the 
1st  Brigade  constitute  the  final  trials:  one  with  only  the  fixed  users 
reporting,  and  one  with  all  units  reporting.  In  preliminary  trials,  the 
effect  of  dropping  those  reporters  which  are  in  another  brigade  was  seen  to 
be  minimal,  except  for  a few  front-line  users  which  had  severe  intervisibility 
problems.  The  errors  for  the  two  primary  runs  are  summarized  in  Table  1 . 

Several  points  deserve  to  be  mentioned  with  regard  to  these  results.  It 
will  be  noted  that  for  30  of  the  44  units,  the  average  error  is  less  when  all 
units  serve  as  reporters.  Furthermore,  in  10  of  the  remaining  14,  the  higher 
average  using  all  units  seems  to  be  due  to  poor  estimates  where  otherwise 
(with  only  the  fixed  units  reporting)  there  would  be  no  estimates.  For  most 
of  the  units  with  maximum  error  above  20m,  those  maxima  occurred  with  less 
than  2 fixed  locators  (i.e.,  where  no  estimate  would  be  available  using  only 
fixed  locators).  The  greatest  error  occurred  for  unit  28,  and  was  108m  j it 


(as  well  as  the  two  unsuccessful  location-operations  for  this  unit)  involved 
21  locators,  including  2 fixed  locators.  (It  is  interesting  to  note  that 
relatively  large  errors  will  occasionally  occur  even  with  a system  " in 
control."  ) Using  only  the  fixed  locators,  that  unit  never  got  lost  and 
never  had  an  error  greater  than  18m  . This  anomalous  behavior  suggests  that 
perhaps  one  should  use  only  the  ranges  from  fixed  units,  if  there  are  enough, 
and  should  use  at  most  a small  number  of  ranges  from  moving  units  in  any  case. 
Also  worthy  of  note  is  the  fact  that  apart  from  this  unit,  every  other 
maximum  error  greater  than  31m  occurred  with  4 or  fewer  locators. 


As  mentioned  earlier,  it  became  obvious  very  quickly  that  low-flying 
aircraft  could  not  be  tracked  from  the  ground.  To  study  the  effects  on 
accuracy  of  such  variables  as  number  of  locators  and  height,  a separate 
study  was  done  on  aircraft.  20  reporters  were  randomly  positioned  on  a 
10  X 10km  piece  of  terrain,  and  kept  fixed.  (See  Fig.  1.)  Ten( simulated) 
aircraft  were  started  at  the  southwest  corner,  at  altituted  of  100,  250, 

400,  550,  700,  1000,  1500,2000,  2500,  and  3500m.  Each  flew  a level  course, 
with  constant  x-velocity;  the  course  reflected  several  times  off  both  north 
and  south  boundaries,  terminating  at  the  eastern  boundary,  systematically 
covering  the  entire  terrain.  (Shown  on  Fig.  1.)  Each  aircraft  went  through 
either  1230  or  2790  location-operations.  The  average  error  and  average 
squared  error  were  tabulated,  both  in  the  x,y-directions  and  in  z,  by 
number  of  locators,  for  each  aircraft  (i9e.,  each  height).  These  tabula- 
tions were  then  studied,  and  regression  techniques  applied. 

The  first  run  counted  locators  for  which  the  range  was  less  than 
8 X height,  12  X height,  16  X height,  etc.,  up  to  80  X he ight, stopping  as • soon 
as  4 or  more  locators  were  found.  It  then  used  LSS  (3d)  . Table  2 shows 
the  numbers  of  trials  for  which  the  algorithm  was  and  was  not  successful  by 
aircraft  height  and  number  of  locators.  There  are  few  cases  where  an  air- 
craft at  700m  or  higher  was  not  located,  but  things  are  not  too  good  below 
that  height.  The  error  statistics  bore  this  out;  the  maximum  z-errors  for 
the  10  aircraft  were  1880,  192,  97,  83,  39,  22,  14,  15,  16,  19m  (in  order  of 
increasing  height),  while  the  maximum  x,y-errors  were  26339,  2373,  337,600, 

72,  36,  29,  23,  27,  18m  respectively.  Since  3d  obviously  does  not  work  with 
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low-flying  aircraft,  an  examination  was  made  of  the  dependence  of  average 
error  on  number  of  locators  and  height,  for  aircraft  heights  of  550m  and 
above.  Regression  analyses  showed  that  the  average  x,y-error  could  be 
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represented  adequately  by  the  expression  E = 16  - ,9n  - — - — + ,017n  , 

x,y  ax  t • 

— 2 

and  average  z-error  by  E ~ 31  - 3n  + .08  n , where  n is  the  number  of 

z 

locators  used.  (Note  that  these  expressions  were  developed  for  aircraft 
higher  than  about  500m.)  The  values  given  by  these  expressions  agree  with 
the  values  obtained  in  the  experiment  within  an  average  deviation  of  0o5m 
for  x,y-error,  and  2m  for  z-error.  (The  experiment  values  for  x,y-error  ranged 
from  3 to  10m;  for  z-error,  from  4 to  29m  .) 


Finally,  a similar  run  was  done  using  slant  range  reduction,  with  air- 
craft at  100,  250,  400,  and  550m  altitudes,  using  p.r.'s,  and  only  those 
reporters  with  elevation  angles  < 20".  The  average  x,y-errors  were  3.2, 

3.3,  3.4,  and  3.7m  respectively;  the  maxima  were  15,  15,  15,  and  17m 
respectively.  This  suggests  that  one  should  use  slant  range  reduction  up  to 
at  least  600m  9 and  use  LSS  3-d  or  slant  range  reduction  above  that  value 
depending  on  whether  there  are  many  good  ranges  at  large  elevation  angles 
or  many  at  small  elevation  angles,  and  on  the  quality  of  altitude  information 
available  at  different  heights.  (It  is  understood  that  p.r.’s  of  the 
accuracy  contemplated  herein  are  limited  in  altitute,  but  will  cover  at  least 
0 to  600m  .) 


Table  1 


Error  surrctmrlos,  Boston  runs. 


All  reporting 

Fixed  units  reporting 

Unit  no. 

Errors  (m) 
Average  Max. 

No . of  times 
not  found 

Errors  (m) 
Avera3e  Max. 

No.  of  times 
not  found 

7 

3.8 

16.6 

49 

5-0 

21.7 

8 

5-9 

39.8 

7-5 

33-4 

105 

9 

8.6 

62.1 

7-0 

39-6 

77 

10 

3.7 

12.0 

• 

5-9 

20.2 

11 

3-4 

9.3 

20 

5.1- 

22.2 

12 

4.2 

16.8 

4.7 

14.2 

13  • 

3-7 

11.2 

104 

5-7 

21.1 

104 

14 

3-5 

10.0 

6.8 

32.4 

15 

3-7 

18.7  • 

4.8 

18.5 

16 

3-7 

23o 

1 

4.8 

16.6 

17 

5.1 

27. c 

• * 

5-1 

25.4 

18 

4-9 

20.4 

5-3 

29.5 

19 

5-0 

18.0 

5-1 

23-3 

20 

3.9 

15.4 

5-8 

30-7 

24 

3.6 

24.1 

4-3 

25.8 

25 

4.3 

13-3 

5-0 

16.3 

73 

26 

2.9 

10.0 

4.1 

11.8 

27 

3-4 

10.1 

4.6 

21.9 

28 

6.3 

108.1 

2 

4.5 

17.7 

29 

3-4 

10-9 

3.4 

15.0 

30 

5-9 

51.6 

4.3 

12.9 

79 

31 

7-7 

79-9 

8 

4.2 

12-3 

61 

32 

3-7 

14.1 

5-1 

33-9 

33 

3-3 

12.  2 

3.6 

12.5 

34 

3-7 

12-5 

3-5  ' 

20.4 

35 

3-2 

11.4 

1 

3-4 

10.5 

36 

3-6 

ll.l 

1 

3-1 

12.8 

37 

3.6 

11.6 

3.3 

13-4 

38  . 

2.7 

10.1 

3.4 

10-7 

42 

3-9 

12.8 

3- .4 

8.5 

102 

43 

3.6 

10.9 

67 

4-5 

16.8 

67 

44 

5-7 

54.8 

3 

4.3 

12.9 

57 

45 

4-9 

23.6 

4.7 

12.5 

69 

46 

3-5 

13.8 

3.8' 

11-7 

67 

47 

9-4 

101.3 

4.0 

11.9 

55 

48 

4.8 

. 27.7 

3 

4.4 

14.2 

77 

49 

3-9 

.12.8 

2 

4.2 

13.4 

2 

50 

8.3 

88.3 

5.9 

26.7 

30. 

51 

5.0 

57-5 

4.4 

19.2 

20 

52. 

3-7 

14.4 

. 5-2 

17-5 

53 

3-8 

13-5 

4.7 

23.1 

77 

54 

3-8 

13-1 

4.6 

25.9 

105 

55 

3-9 

17.2 

5.4 

25.5 

56 

4.2 

16.4 

5.7 

22.0 

All 

4-5 

108.1 

4.7 

39-6 

• 

Table.  2 


Lost  and  found  summary  aircraft  study 
Aircraft  height 


Number  of 
Locators 

100 

250 

400 

550 

700 

1000 

1500 

2000 

2500 

3500 

3 

4/0 

4 

136/640 

84/396 

44/175 

18/80  - 

7/52 

5 

40/255 

48/233 

28/189 

18/92 

3/20 

6 

18/98 

18/163 

8/131- 

9/103 

0/11 

7 

7/24 

14/98 

22/192 

16/128 

4/32 

0/1 

8 

1/5 

6/57 

4/122 

2/119 

0794 

0/2 

9 

0/2 

6/83 

3/189 

5/159 

1/57 

0/1 

10 

1/22 

2/95 

1/115 

1/70 

0/3 

11 

0/1 

0/16 

0/85 

1/98 

0/0 

12 

0/10 

0/94 

1/130 

1/28 

13 

0/90 

0/44 

0/13 

14 

0/70 

0/65 

0/36 

15 

0/16 

2/228 

0/184 

16 

0/10 

0/133 

0/154 

17 

0/44 

0/212 

18 

0/33 

0/231 

19 

0/232 

0/67 

20 

0/132 

1/1162 

2/1228 

0/1230 

0/1230 

Total 

Lost 

206 

177 

111 

69 

19 

1 

1 

2 

0 

0 

Fig.  1 

Locators  1 positions  and  aircraft  path.  Aircraft  error  study. 


4.  Conclusions 


The  following  represent  conclusions  with,  varying  degrees  of 
certainty;  some  only  say  what  might  be  true^  or  what  should  be  checked 
out.  They  are  organized  to  parallel  the  earlier  list  of  questions. 

(1)  The  system  will  work  under  reasonable  conditions,  such  as 
outlined  below. 

(2)  The  system  will  generally  locate  ground  units  within  reasonable 
error  bounds  whenever  such  units  have  at  least  4 fixed  locators 
or  well-located  moving  locators  (which  may  be  airplanes,  as  long 
as  they  are  not  too  nearly  overhead) . 

(3)  x,  y-errors  will  be  the  same  order  of  magnitude  as  range  errors, 
for  about  4-6  good  locators;  the  accuracy  (error)  should  improve  as 
l//n-2  where  n is  the  number  of  locators,  up  to  the  point  where 
the  computer  word  length  begins  to  limit  the  algorithm !s  accuracy. 
The  z-errors  will  be  comparable  for  genuine  3-d  situations;  they 
will  depend  solely  on  the  accuracy  of  the  independent  pressure 
reference  when  that  is  used. 

(4)  Use  LSS;  to  find  a starting  value  when  extrapolation  from  previ- 
ous estimates  is  not  to  be  trusted,  use  LSL.  Use  an  outlier^ 
rejection  scheme  in  each,  in  addition  to  a preliminary  screening 
of  the  incoming  ranges  to  eliminate  gross  (10%)  errors.  Use  the 
2-dimensional  version  on  ground  units,  with  slant^range  reduction 
if  the  terrain  is  not  smooth  (unless  all  but  1 or  2 of  the  locators 
are  aircraft  at  high  elevation  angles,  say  60°  or  more, in  which 
case  3-d  treatment  is  required) ; also  use  2-d  with  slant-range 
reduction  for  aircraft  with  altitudes  less  than  about  1000m, 
possibly  in  conjunction  with  3d  for  aircraft  between  (say)  600 

and  1000m.  Use  3d  for  aircraft  over  1000m.  (Note;  These  numbers 
depend  on  the  spread  of  ground  units ; the  values  given  are  appro- 
priate for  the  Boston  deployment.) 


(5) 


(6) 

(7) 

(8) 


(9) 


(10) 


(11) 


(12) 


Fixed  units  should  constitute  a large  fraction  (at  least  1/2)  of  the 
set  of  reporters.  The  system  will  work  with  fewer.  But  the  possibil- 
ity exists  that  the  set  of  users'  estimated  positions  will  "break  its 
moorings"  and  wander  away,  as  one  entity,  from  the  true  position 
It  is  not  necessary  that  the  same  units  be  fixed  throughout  the  period 
Of  time  that  the  system  is  in  use.  (In  which  case  it  follows  that 
the  positions  in  which  they  stop  will  not  be  "known",  only  estimated.) 

This  question,  concerning  the  number  of  reporters  needed,  is  answered 
under  (3)  and  (4)  above. 

It  is  understood  that  maps  will  not  be  generally  available,  but  this 
is  fortunate  in  a way,  since  their  use  extracts  a very  heavy  penalty 
m computation  time.  Altimeters  (or  more  accurately,  "pressure  ref- 
erences", or  p.r.'s)  will  be  necessary  when  aircraft  at  low  elevation 
angles  are  working  with  ground  units,  either  as  locatee  or  as  locators 
and  will  be  helpful  whenever  the  terrain  is  rough  enough  that  altitude 
differences  among  ground  units  are  as  large  as  the  error  a of  the  p.r 
accuracy  .of  14m  (a)  in  p.r.  values  (that  is,  a = 20m  for  the 
difference  between  the  heights  of  two  units)  will  enable  the  system 
to  maintain  position  errors  comparable  to  range-measurement  errors. 
Better  p.r. 's  will  of  course  permit  more  accurate  position  estimates. 

eporting  once  per  30sec.  seems  reasonable,  except:  a)  Moving 

locators  should  report  more  frequently,  especially  aircraft;  and 
b)  if  the  system  is  to  be  used  for  aircraft  guidance,  close  support 

or  fire  control,  then  all  aircraft  should  report  more  often,  say  once 
every  second  or  two. 

To  use  genuine  3d,  a minimum  of  3 locators  is  required;  only  one  need 
be  an  aircraft.  However,  at  least  5 (including  2 aircraft)  is  advisable 
so  that  bad  measurements  may  be  detected.  Note  that  aircraft  at  low 
elevation  angles  count  as  ground  units  here, 
and  (13)  - See  answer  to  (4). 


(14)  Use  as  many  locators  as  possible,  but  properly  weighted. 

(15)  Weight  fixed,  well-known  locators  considerably  higher  than  others. 

(See  discussion  of  weights  earlier  in  this  section.) 

(16)  Continual  relocation  is  not  advisable.  Periodic  relocation  may  be 
worthwhile;  compare  with  (19).  below. 

(17)  Do  not  use  only  fixed  units  as  reporters,  except  in  situations  with 
unusually  good  intervisibility.  (There  are  too  many  times  when  a 
unit  won’t  have  enough  fixed  locators.) 

(18)  There  are  indications  that  relatively  large  errors  will  occur  when 

a locator  begins  to  move,  if  the  system  continues  to  treat  his  posi- 
tion as  fixed.  Also,  whatever  time  it  takes  to  establish  that  a unit 
has  stopped  becomes  a delay  in  "fixing"  that  unit’s  position  so  that 
he  can  become  a good  locator.  Therefore  a definite  penalty  is  attached 
to  not  having  a moving/stopped  indicator.  The  size  of  this  penalty 
can  only  be  established  after  the  appropriate  filtering  techniques  are 
implemented.  Thus  a complete  answer  awaits  a development  study. 

(19)  Indications  are  that  it  is  possible  to  resurvey  the  set  of  fixed  units. 

In  the  course  of  system  development,  attention  should  be  paid  to  devel- 
oping, implementing,  and  testing  a multivariable  least-squares  approach, 
allowing  the  adjustment  of  all  but  a few  of  the  fixed  reporters’  posi- 
tion estimates  simultaneously.  This  technique  would  be  applied  peri- 
odically during  long  battles,  when  the  set  of  currently-fixed  reporters 
becomes  considerably  different  from  the  set  of  initially-fixed  reporters. 
(This  technique  might  be  either  a supplement  to,  or  a replacement  for, 
the  periodic  use  of  relocation;  see  (16)  above.) 
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ABSTRACT 


The  analysis  of  error  propagation  for  position  location  in  a one- 
dimensional geometry,  assuming  Gaussian  error  distributions  and  maximum- 
likelihood  estimation  (i.e.,  location  by  minimizing  a particular  quadratic 
penalty  function),  is  mathematically  tractable.  Although  this  scenario  is 
unrealistic,  "solving"  it  in  closed  form  may  yield  insights  applicable  to 
more  general  situations.  Working  Paper  No.  4 developed  formulas  for  the 
law  of  error  propagation  through  a single  location-operation  in  such  a 
situation.  The  present  note  proceeds  to  treat  a series  of  location-opera- 
tions for  the  (simplest)  case  where  only  two  units  are  involved,  and  shows 
that  in  this  case  the  precision  of  location  stabilizes  over  time.  It  is  • 
planned  next  to  seek  to  extend  this  encouraging  result  to  more  than  two 
units. 


Note : Project  working  papers  are  informal  documents  prepared  to  facilitate 

they  may  contain  tentative  or  relatively 


discussion  and  communication; 
unchecked  material. 


0.  INTRODUCTION 


Following  the  terminology  of  other  papers  in  this  series , we  shall 
use  the  term  location-operation  to  refer  to  the  process  of  imputing  a 
position  to  one  unit  (the  locatee)  , on  the  basis  of  its  measured  distances 
from  a set  of  other  units  (the  observers)  whose  estimated  positions  are  at 
hand.  Both  the  distance  measurements ; and  the  estimates  of  the  observers' 
locations,  are  subject  to  error.  The  purpose  of  Project  "WHERE"  can  be 
roughly  described,  from  a technical  viewpoint,  as  that  of  determining  (for 
a variety  of  scenarios)  the  propagation  of  location-error  over  a series  of 
such  operations,  in  which  the  role  of  "locatee"  varies  cyclically  among  a 
large  set  of  units,  many  of  which  may  be  in  motion  during  part  or  all  of 
the  series.  As  a principal  tool  for  this  purpose,  project  staff  are 
developing  the  various  components  of  a computerized  simulation  ("WHERESIM") . 

While  such  a simulation  is  needed  for  the  examination  of  error-propaga- 
tion in  "general"  or  "unpatterned"  situations,  it  is  plausible  that  valuable 
insights  can  be  obtained  from  the  study  of  special  situations  which  are 
simple  enough  to  admit  direct  mathematical  analysis.  Working  Paper  No.  4^^ 
initiated  one  such  study,  whose  key  simplifying  assumptions  are  those  of 

(a)  one -dimensional  geometry,  and 

(b)  Gaussian  error  distributions. 

Specifically,  WP4  derived  the  error-propagation  formula  under  those  assump- 
tions for  a single  location-operation,  leaving  for  a next  phase  of  work  the 
recursive  application  of  this  formula  to  determine  error  build-up  through  a 
sequence  of  location-operations,  with  various  units  serving  successively  as 
locatee.  The  present  paper  carries  out  this  further  analysis  for  an  especially 
simple  situation,  that  in  which  only  two  units  are  involved.  The  units  thus 
alternate  between  the  roles  of  locatee  and  observer. 


Project  "WHERE"  Working  Paper  No.  4,  Error  Propagation  for  Position 

Location  in  One  Dimension,  M.H.  Pearl,  December,  1971.  To  be  referred  to 
in  what  follows  as  "WP4." 


The  analysis  in  WP4  was  divided  into  two  cases,  according  as  the 
single  location-operation  considered  there  did  or  did  not  make  use  of  a 
prior  estimate  of  the  locatee's  position.  From  the  results  in  WP4  (see 
its  last  page)  and  a variety  of  additional  evidence  developed  in  other 
parts  of  Project  "WHERE"  , it  now  seems  clear  that  such  prior  estimates 
should  indeed  be  employed,  and  so  their  use  is  assumed  in  what  follows. 
Specifically,  each  such  prior  estimate  is  obtained  from  (i)  the  unit's 
estimated  position  at  the  end  of  the  last  previous  location-operation  and 
( ii)  its  velocity  of  motion(assumed  known  exactly)  in  the  interim. 

With  the  decision  between  "prior-using"  and  "prior-ignoring"  thus 
made  in  favor  of  the  former,  the  analysis  of  the  present  paper  is  split 

(2) 

into  two  cases  along  a different  line  of  cleavage.  In  Working  Paper  No. 2 , 

it  was  suggested  that  the  measurements  and  calculations  taking  place  during 
a location-operation  might  be  exploited  not  only  to  locate  the  locatee,  but 
also  to  revise  the  currently  assumed  positions  of  the  observers,  so  that 
one  would  in  effect  have  a "location-and-revision"  operation.  Accordingly, 
both  "with  revision"  and  "without  revision"  scenarios  will  be  analyzed 
below.  The  most  significant  finding  is  that-  (in  both  cases)  the  precision 
of  location  stabilizes  over  time.  (The  natural  next  step  is  to  investigate 
whether  this  encouraging  result  combines  to  hold  for  more  than  tv?o  units.) 

Although  familiarity  with  WP4  is  not  essential,  it  will  be  quite  helpful 
in  reading  the  present  text,  which  is  organized  as  follows.  Section  1, 
appearing  next,  describes  the  underlying  mathematical  models  involved  and 
sets  up  most  of  the  necessary  notation.  Section  2 presents  the  error-propaga- 
tion analysis  for  a sequence  of  "location-with-revision"  operations,  while 
Section  3 does  the  same  for  the  "without  revision"  case. 


(2) 


Project 


"WHERE" 


Working  Paper  No.  2,  An  Iterative  Approach  to  the 


Location-Operation , A.J.  Goldman,  October,  1971. 


. 


1.  MATHEMATICAL  FORMULATION 


The  analysis  deals  with  2 units  located  on  a line  (i.e.,  the  geometry 
is  one-dimensional),,  The  relevant  time  period  begins  at  time  t * 0 , and 
the  successive  location-operations  occur  at  times  t.  (j  - 1 ,2 , ...)  where 


0 = t < t.  < t0  < . . . . 

o 1 2 

For  i = 1,2  and  j = 0,1,2,...,  we  set 


so  that 


U. , = true  position  of  unit  i at  time  t. 

ij  J 


D. 

J 


true  (signed)  distance  between 
the  units  at  time  t 


(1.1) 


The  1 s are  unknown  to  the  location-system,  which  jls  however  "given” 

at  time  t.  quantities. 


u..  — prior  estimate  of  U. . at  time  t.  , 
ij  iJ  J 


to  be  discussed  further  below.  The 
tion-operation  is  based  consists  of 


information  on  which  the 


u2 j 3 and 


j-th  loca- 


d.  = measured  value  of  D . 

J J 

This  measurement  is  assumed  unbiased,  i.e.  E(D.)=d,  , and  also  normally 

J J 

distributed,  so  that  the  random  measurement  error 


S . *=  d . - D . 

J . J J 


(1.2) 


■ 


is  Gaussian  and  satisfies  (for  some  a > 0 ) 


E(S  ) = 0 , E(52)  = a2  . 


(1.3) 


Note  that  our  notation  involves  a and  not  O.  ; i.e.,  we  assume  no 

J 

systematic  change  over  time  in  the  quality  of  the  distance  measurements. 


The  initial  position-estimates,  u^q  and  u^q  , are  assumed  to  be 

unbiased.  i„e.  E(u.  ) = U.  , and  also  to  be  normally  distributed;  in 
3 10  10  3 3 

terms  of  the  vector  u defined  by 

o 


T 

u 

o 


(u10>  u20^  ’ 


we  set 


Q = var(u  ) 
o N o 


var(u10>  COV^U10,U20^ 


COV(U10'U2(P  var(u2o^ 


(1.4) 


the  matrix  Qq  (describing  the  quality  of  our  start ing  position-data) , as 
well  as  a (which  describes  the  quality  of  distance  measurements),  are  the 
basic  inputs  to  our  analysis  of  error  propagation. 


It  is  assumed  that  between  t..  , and  t.  , unit  i moves  with  known 

j-i  j 

velocity  v„  ; it  therefore  covers  a known  distance 


A.  . 


vij(tj 


(1.5) 


Thus  the  prior  estimates  of  the  two  units'  positions,  just  before  the  first 
location-operation,  are 


u 


il 


u . + A. 

10  . 10 


(i  = 1,2)  ; 


it  follows  from  (1.4)  that 


Qq  = var(Ul)  , 


(1.6) 


where  we  have  used  for  j = 1 the  notation 


u.  = (u  _,u0 . ) . 
J Ij  2j 


Let  us  define  the  vector  by 


U 


<W  • 


The  calculations  of  the  j-th  location-operation,  in  our  present  model , 
yield 


U.  — maximum- likelihood  estimate  of  U. 
J J 


based  on  the  data  u.  and  d.  ; 

j j ’ 


this  estimate  is  known  (cf . WP4)  to  be  unbiased,  and  we  put 


Q.  = var (U. ) 
J N J 


(1.7) 


It  is  at  this  point  that  the  "with  revision"  and  "without  revision" 

-k 

cases  diverge.  In  the  first  case,  U.  .is  accepted  as  the  vector  of  loca- 

1 * 

tion  estimates  after  the  j-th  location-operation;  thus  represent  the 

quality  of  our  position- locations  at  time  j , so  that  our  aim  is  to 

•k 

determine  these  matrices  (and  in  particular,  their  behavior  as  j 

increases).  Moreover,  in  this  case  the  "next"  set  of  prior  estimates  is 
given  by 


u .,i  = U + A.  , 

j+i  j i ’ 


(1.8) 


where  we  have  defined  A.  by 

J 


■ 


In  the  second  ("without  revision") 
after  the  j-th  operation,  denoted 
is  taken  as 


case,  the  vector  of 

kk 

U.  with  components 


location  estimates 


Q 


kk 

ij 


and  U 


2j 


kk 


kk 


Ulj  Ulj  J U2j 


kk  kk 

Ulj  = Ulj  ’ U2j 


k 


if  j odd  (unit  1 


if  j even  (unit  2 


the  locatee)  , 

(i. 

the  locatee)  , 


so  that  the  prior  estimate  of  the  observer’s  position  indeed  goes  unrevised. 
The  "updating"  formula  (1.8)  is  of  course  replaced  by 


And  now  the  matrices 


kk 


u.-t  = U.  + A. 
J+l  J J 


kk  _ kk % 

Q.  = var(U.  ) 
J J 


(1.10) 


(1.11) 


are  the  matrices  which  are  to  be-  investigated,  and  whose  limiting  behavior 
(for  large  j ) is  to  be  investigated. 

The  analyses  of  the  two  sequences  [Q^  J and  (Q^.  J of  matrices,  are 
the  topics  of  Section  2 and  3 respectively. 


c? 


2.  ANALYSIS  ASSUMING  REVISION 


It  is  convenient  to  define  the  random  error  vectors 


e.  = u.  - U. 

J J J 


(2.1) 


Then  the  likelihood  function  for  the  first  location-operation  is 


L = K 


exp{-%(ej[  Qg1  ex  + a'2  &2  ) } 


where  K is  a positive  constant,  so  that  maximizing  L is  equivalent  to 
minimizing 


"k  T — 1 — 9 9 

L = %(«£  Q0  \ + a z Sp  . 


(2.2) 


From  (1.1)  and  (1.2)  , 


b,  = d.  + [1,  -1 ]U-  , 


and  from  this,  (2.1),  and  (2.2),  it  follows  that  the  maximum- likelihood 
* 

estimate  is  the  value  of  which  minimizes 


L*  = %{(ux  - Up1  Q"1^  - Ux) 

+ ^‘2(d1+  [1,-1]  u1)T(d1+  [1,-1]  Uj_)  } . 


(2.2a  ) 


We  next  apply  the  (straightforwardly  verified)  formula 


3 [(AX  + B)T  S(AX  + B)]/3X  = 2AXS(AX  + B) 


(2.3) 


where  X is  a column  vector  and  S is  a symmetric  matrix.  Taking  A - I 
(the  2 X 2 identity  matrix),  B = u^  and  S = Q“^  , 


we  obtain 


' 


(2.4) 


3{VVT  Qo1(urUl)}/5Ul  = 2Q01(Ul"ul)  5 

taking  A 55  [1,-1]  , B = d^,  and  S = I ; we  obtain 

3{(dx+  [l,-l]U1)T(d1+[l,-l]U1)}/3U1  = 2[l,-l]T(d]+[l,-l]U 

= 2{  [1,-1  ]T  dj+  NUX} 

where  we  have  put 


N = [1,-lf  [1,-1] 


1 

-1 


Thus  the  equation  dL  /dU^  **  0 , to  be  satisfied  by  **  , 

Qo1(urui>  + a’2t1^-1]  di  + a'2  nux  = 0 ; 

setting 

H *=  a”2  Q N , 

o o’ 

we  see  that  this  is  equivalent  to 


(I  + Hq)U*  = ux  + a'2  Qq[1,-1]  dx  , 

so  that 

U*  = (I  + Ho)'1  (ux  + a'2  Qo[l,-l]  dx)  . 
Since  cov(u^d^)  = 0 y we  have  by  (1.6)  and  (1.8) 

var(.uX  + CT"2  dp 

- Q0  + o'4  Q0[l,-l]a2  [1,-1 f Qc 


) 

, (2.5) 

(2.6) 

reads 


(2.7) 


(2.8) 


(2.4) 


^VV1  Q^CvVJ/aUjL  “ ; 


taking  A a [1,-1]  , B - d^,  and  S = I , we 


obtain 


3f(d1+  [l,-l]U1)T(d1+[l,-l]U1)}/SU1  = 2[l,-l]T(d1+[l,-l]U] 

= 2{  [1,-1  ]T  dj+  nux5 


where  we  have  put 


n = [i,-ir  [i,-i]  « 


i -i 

-i  i 


"k  . ~k 

Thus  the  equation  dL  /dU^  = 0 } to  be  satisfied  by  *=  , 


Q"1(U1-U1>  + a'2[l,-l]  dx  + a"2  HUX  - 0 ; 


setting 


H = a“2  Q N . 

o o 7 


we  see  that  this  is  equivalent  to 


(I  + Hq)U*  = ux  + a*2  Qq[1,-1]  dL  , 

so  that 

U*  = (I  + Ho)_1  (ux  + a-2  Q0[l,-1]  dp  . 

Since  cov(u^^d^)  = 0 } we  have  by  (1.6)  and  (1.8) 

var(ux  + a'2  Qo[l,-l]  d^ 

“ Qo  + CT’4  ^o[1'"1]a2  [1^1]T  Qo 


) 

, (2.5) 

(2.6) 

reads 


(2.7) 


(2.8) 


(2.9) 


= Q + cr 
o 


-2 


Q N Q = (I  + H )Q 
o o o o 


Since  (2.8)  yields 


Q1  = var(U1  ) 

= (I  + Ho)_1  var(Ul  + a'2  Qo[l,-l]d1)(I  + Ho)"T  , 
it  follows  from  (2.9)  that 

•k  -T 

Q1  - VI  + V > 

which  is  equivalent  to 

< = (I+H0)_1Q0  . 


We  turn  now  to  the  second  location-operation.  By  (1.8), 


"k  k 

var(u2)  = var(U1)  = Qx  , 
k 

i.e.  plays  the  same  role  now  that  was  played  by  previously, 

if  by  analogy  with  (2.7)  we  put 

Hx  = a 2 N , 

then  by  analogy  with  (2.10)  and  (2.11)  we  will  obtain 


Q2  = V1  + 


k -T  * -1 

H,)  «*  (X  + HO 


(2.10) 


(2.11) 


Thus , 


. 


to 


is  given  by  the  two 


More  generally,  the  transition  from 
equations 


Q 


kk 

j+1 


* 

H. 

J 


-2  * 

* qj 


N 


• k k , * -T 

Q.  = Q.(I  + H.) 
J+l  J J 


* -1  k 

(I  + H.)  Q. 
J J 


(2.12) 

(2.13) 


k 

Q - Q 

o o 


with  initial  conditions  0 =0  and 

Using  these  results,  we  will  prove  by  induction  on  j 


* 

H = H 
o o 


that 


Q*  = (I  + jH  )'1  Q = Q (I  + jH  )"T 

J v o o o o 


(2.14) 


By  (2.11),  it  is  true  for  j = 1 . Suppose  it  is  true  for  some  particular 
value  of  j ; we  will  prove  it  correct  for  j + 1 . First,  by  (2.12)  and 
(2.14), 

JL  A i i i 

H.  = a"z(I  + jH  ) Q N = (I  + jH  ) H = H (I  + jH  ) . (2.15) 

J v o o o o o o 

If  we  set 

Z.  = (I  + H*  )_1  , 

J J 

then 

I 


so  that 

(I  + jH  ) « Z . (I  + jH  ) + Z ,H  - Z . [I  + ( j+l)H  ] 

O J O JO  J O 

and  thus 

= (I  + jHo)[I  + (j+DHj-1  = [I  + (j+l)H0]_1(l  + jHo); 


= Z . (X  + H*  ) , 

J J- 

= Z.  + Z.H  (I  + jH  ) 
J jo  oy 


Z. 

J 


it  follows  from  (2.13  - 15)  that 


* 

Vi 


* 

Z.Q. 

J J 

[I+(j+l)Ho]_1[I  + jHQ][l  + jHo]-1  Qo 
[X+(j+l)Ho]-1  Qq  , 


so  that  (2.14)  also  applies  for  j-fl  . The  induction  proof  is  complete. 


Note  that  with  (2.6),  (2.7)  and  (2.14),  we  have  succeeded  in 

k 

determining  the  matrices  explicitly  in  terms  of  the  basic  data 

and  ex  . We  now  prove  that  the  limit 


= lim 


CO 


* 

Qj 


(2.16) 


exists,  i.e.  that  the  precision  of  position-location  stabilizes  over  time, 

k 

and  moreover  we  shall  evaluate  explicitly. 


Let 


then  by  (2.6-7) 


I + jH  * 
J o 


1 + cr  2j(qx-q)  o 2j  (q  - 


-2. 
<7  J 


-2 . / N 
o J(q  ~ q2> 


(a  /j)  + (qx-q) 

(q  ~ q2) 


-2 

1 + O j(q0  - q) 


q - qi 


(o  /j)  + (q2-q) 


(2.17) 


■ 


The  coefficient  of  a is  a matrix  with  determinant 


(a4/j2)  + (a2/j)(q1  + q2“2ci)  = (tf2/j)[(a2/j)  + (qx+  q2~  2q)] 


and  thus  with  inverse  matrix 


(a2/j)  1 [(a2/j)  + (qx+  q2~  2q) ] 1 


(a2/j)  + (q0“  q) 


qx-  q 


q2“  q (<?  /j)  + (q^  q) 


It  follows  from  (2.17)  that 


(I  + jHQ)  1 = [(cf2/j)  + (q^  + q2-  2q)]  1 


(c  /j)  + (q2  - q)  qx-q 


q2-  q 


(a2/j)  + (qrq) 


In  particular. 


(2.18) 


lim  . (1  + jH  )_1  = H* 

J -»■  CD  N O 


(2.19) 


exists  and  is  given  by 


H s (qx  + q2  ~ 2q)  1 


q9  - q qi  - q 


q9  - q 


«i  - q 


(2.20) 


. 


thus  by  (2.14)  the  limit  exists  and  is  given  by 


= (q1+q2‘2<i>'1  (q1q2‘ci2)  J > (2-21) 

where  J is  a 2X2  matrix  of  ones.  It  is  somewhat  surprising  that  the 
limit  is  independent  of  a 

A natural  question  is  whether  the  position-locations  grow  increasingly 

* * 

precise  over  time,  i.e.  whether  the  diagonal  entries  of  = var(tL  ) are 

decreasing  functions  of  j .By  (2.14)  and  (2.18),  the  (1,1)  entry  of 
* 


var(U^j)  = q1[(a2/j)  + q^Uj^"  q2)  3/ [ Cct2/ j ) + (q1+q2"2q)]  • (2.22) 

Now  it  is  readily  verified  that  the  function 

f (x)  s»  (x+a)/(x+b)  = 1 - (b-a)/(x+b) 

is  increasing  or  decreasing  according  as  (b-a)  is  positive  or  negative. 
Applying  this  with 

2 -1  2 

x = a /j  , a = qx  (q^  - q ) , b 88  qx  + q2  “ 2q  , 


ql  + q2  “ 2q  ~ ql1^qlq2  " q2^  ^ • (2.23) 

2 

which  indeed  holds  since  it  is  equivalent  to  (q.  - q)  > 0 , Analogously, 

* 

var(U2j)  is  a non- increasing  function  of  j 


' I :m 


3.  ANALYSIS  FOR  LOCATION  WITHOUT  REVISION 


We  now  turn  to  the  matrices 


var(U.  ) 
J 


(3.1) 


corresponding  to  the  location-without-revision  scenario.  The  analysis  is 
more  complicated  here,  since  the  two  units  no  longer  figure  symmetrically 
in  each  location-operation  (recall  that  unit  1 is  the  locatee  in  odd-numbered 
operations,  unit  2 the  locatee  in  even-numbered  ones). 

It  is  convenient  to  repeat  from  Section  2 the  definitions 


column  k ) position,  and  all  other  entries  zero.  Thus,  for  example,  the 
consequence 


N = [1,-1]T  [1,-1] 


1 -1 


(3.3) 


-1  1 


H 


o 


(3.3) 


of  (1.9)  can  be  written 


■ • ! 


From  this  and  (3.1)  it  follows  that 


irk 


Qx  = En  var(Ux)  En  + E22  va^u^E^ 

+ EX1  cov(U1 , ^1)E22  + E22  cov(ui;  Ux  ) E1X 


(3.5) 


By  (2.8)  and  (2.11) 


cov(U1^u1) 


(I  + H )“1Q 
o o 


- Q,  > 


so  that  (3.5)  yields 


Q1  E11  Q1  E11  + E22  Qo  E2: 


+ EX1  Ql  E22  + E22  Qx  E1X 


= <*1  + E22<V  ^1>E22 


By  use  of  (2.11),  this  reads 


Q**  = (I  + Ho)'1  Qo  + E22[I  - (I  + Ho)'1]Q„  E 


o 22  ' 


(3.6) 


The  second  location-operation  follows  the  same  logic  as  the  first, 
irk 

except  that  plays  the  part  of  Qq  and  that  the  roles  of  the  two  units 

are  interchanged.  Thus,  by  analogy  with  (3.3)  and  (3.6),  we  have 


H. 


irk 


a’2  Q,  N 


irk  irk^  - ] irk  , irk  -1  _ irk 

>2  = (X  + Hj_  ) Qx  + Eutl  - (X  + H1  ) ] Qx  En 


(3.7) 

(3.8) 


tT 


More  generally,  the  appropriate  recursion  equations  are 


kk  - 9 ** 

H.  = cr  Q.  N 
J J 


kk 

Q2j+1 


kk  -1  kk  kk^  -1  kk 

(I  + H2j)  Q2j  + e22[i  - (I  + H2j)  ]Q2.  e22  , 


(3.9) 

(3.10) 


Q2j+2  = (1  + H2J+1)  Q2j+1  + En[X  - (I  + H^)-1]  Q2j+1  Eu;  (3 


these  are  accompanied  by  the  initial  conditions  Q =0  and  H ~ H 

o^ooo 


In  view  of  the  way  in  which  the  analysis  "alternates"  between  the  two 
units,  it  would  be  unreasonable  to  expect  the  type  of  convergence  given  by 
the  existence  of  a single  limiting  matrix 


kk  kk 

Q = 1 im . Q . 

^00  J •>  00  ^ j 


Rather,  one  might  hope  for  "convergence  with  period  2"  , i.e.  for  the 
existence  of  a pair  of  matrices 

Q°  = lim.  Q0.  - , Q6  ~ lim  Q0.  , (3.12) 

\»  j ->  CD  X2j+1  ^ OO  j ->  CD  ^2j  3 

where  the  superscripts  "o"  and  "e"  stand  for  "odd"  and  "even" 
respectively.  1 Furthermore,  in  view  of  the  results  obtained  in  Section  2,  one 
would  not  be  surprised  if  these  two  limiting  matrices  turned  out  to  be 
independent  of  a . 

At  present  we  have  not  succeeded  in  analyzing  this  situation  with  any- 
thing like  the  completeness  achieved  in  Section  2 for  the  "location  with 
revision"  case.  However,  we  have  succeeded  in  proving  that  in  the  present 


F 


$X 


■ 


case,  like  that  of  Section  2,  the  precision  of  position-location  stabilizes 
over  time.  To  make  this  precise,  let 


Qj 


(j) 

(j) 

1 

q 

(j) 

(j) 

then  = q and  q^°^  = q^  (i  = 1,2)  . The  precision  with  which  the  j-th 

location-operation  estimates  the  position  of  the  locatee  is  given  by 


a) 

2 

= var (U*  ) 

if 

j is  even  , 

(j> 

1 

= var(U, .) 

if 

j is  odd 

It  will  be  shown  below  that 


sequences  {q^^}  and  are  convergent,  (3.13) 


thus  proving  the  "stabilization"  mentioned  above. 


The  logic  of  the  proof  will  run  as  follows.  It  will  be  shown  that 


sequences  {q^*^  } and  {qjj^*^  ] are  convergent  . 


(2j) 


(3.14) 


This  is  a property  of  the  sequence  {q  } of  even-index  matrices.  Since  the 

sequence  {Q 2 j-j-i 5 °f  odd-index  matrices  obeys  the  same  type  of  recursion, 
and  since  (3.14)  is  symmetric  as  between  the  two  units,  it  follows  by  analogy 


that 


sequences 


{qi(2J+1)} 


and 


r (2j+l)-> 

Uo  i 


are  convergent . 


(3.15) 


The  desired  result  (3.13)  then  follows  from  (3.14)  and  (3.15). 


Thus  it  suffices  to  prove  (3.14).  Since  the  variance  terms  q 


(2j) 


and 


(2j) 


are  nonnegative,  i.e.  bounded  below  by  0 , it  suffices  to  prove  that 


the  sequences  in  (3.14)  are  non- increasing . These  statements. 


q(2j+2)  < q(2j) 


(2j+2) 


< q 


(2j) 


are  to  be  proved  by  induction  on  j , but  it  suffices  to  establish  the 
prototype  induction  step,  namely 


(2)  . „ 
ql  ^ ql 


(2)  < n 
q2  - q2 


(3.16) 


The  proof  of  the  desired  result  (3.13)  has  now  been  reduced  to  the 
demonstration  of  the  two  statements  (3.16)  . For  proving  them,  it  is 
convenient  to  introduce  the  auxiliary  quantities 


s = (qx-q)  + (q2-q) 


A 


2 

qlq2"q 


(3.17) 


Since  Q is  a variance-covariance  matrix,  its  determinant  A > 0 , and 
it  follows  by  (2.23)  that  S > 0 . The  following  calculations  will  make  use 
of  the  easily-checked  identities 


A - qs  = (q1-q)(q2“ci)  * 

2 

A - = -(q^q) 


(3.18) 

(3.19) 


I 


Direct  calculation  from  (3.6)  yields 


**  ,2  N-1 

Qn  - (a  + S) 


2 2 
a + A a q + A 


2 2 
a q + A (a  + S)q, 


(3.20) 


From  (3.7 ),  we  obtain 


**.-1  r 4 2 N 2,-1 
(I  + Hx  ) 1 - [c  + 2Sa  + (q2«q)  ] 


a4  + (S  + q2-q)a2  + (q2*-q)2  ^(q^^-q) 


ff2(q2“q)  + (q2"q^2 


a +(S+q--  q)cr 


Application  of  (3.8)  and  (3.20)  yields 


= (a2  + S)‘1(a2q1  + A)  , 


so  that  the  first  of  the  desired  statements  (3.16)  is  equivalent  to 


cr  q^  + A < (a  + S)q^ 


and  hence  to  A - q^S  < 0 , which  is  true  by  (3.19).  Next,  substitution 

into  (3.8)  yields 

qo2^  = (cr  + S)  [a4  + 2Sa2  + (q0-q)2]  1 {q9a^  + [A  + 2q_S  ]a4 


+ [2q  S2-  (q9-q)2(q9-2q) ]a2  + A(q9-q)2}  , 


. 


■ 


so  that  the  second  of  the  desired  statements  (3.16)  is  equivalent  to 


q2a6  + [A  + 2q2S]a4  + [2q,,S2-  (q2-q)2(q2-2q)  ]a2  + A(q2-q)2 

< q2^CT^  + ^sa  + Cq2-<i)2 1 i 


which  is  true  since  it  is  in  turn  equivalent  to 

(q2-q)2  [C2  + (q2-q)]2  > 0 . 

This  completes  the  proof  of  (3.13),  so  that  the  existence  of  the 
"variance  terms"  in  (3.12),  i.e.  entries  q°  and  q°  in  Q°  and  entries 
e e e 

q^  and  q^  in  , is  established.  We  have  not  as  yet  succeeded  in  proving 

the  existence  of  the  covariance  limits 


o 

q 


(2j+l) 

->  CD  ^ 


lim. 

J 


(2j) 


->  OO 


(3.21) 


i.e.  the  off-diagonal  terms  in  (3.12),  except  under  the  additional  hypotheses 


qx  > q ; q2  > q * (3.22) 

This  assumption  is  probably  not  very  restrictive  for  the  situations  of 
interest  here;  it  asserts  that  the  errors  in  the  initial  estimates  (of  the 
two  units’  positions),  if  positively  correlated,  have  covariance  less  than 
either  of  their  individual  variances.  Because  A > 0 , (3.22)  will  automat- 
ically be  satisfied  if  q^  = . 


To  prove  that  (3.22)  implies  (3.21),- 
ilization 


we  first  prove  that  it  implies  the 


(3.23) 


■ 


i?fo 


of  itself.  The  proof  is  by  induction  on  j , and  since  (3.23)  is  symmetric 
as  between  the  two  units , it  suffices  to  treat  a prototype  induction  step, 
i.e.  to  show  that 


> q 


(i) 


(3.24) 


follows  from  (3.22).  From  (3.20),  we  see  that  the  first  part  of  (3.24)  is 
an  immediate  consequence  of  the  first  part  of  (3.22).  Also  from  (3.20),  we 
see  that  the  second  part  of  (3.24)  is  equivalent  to 


(a2  + S)q9  > a2q  + A , 


hence  to 


2 2 

o’  (q9-  q)  > A - q9s  - (q9-  q)  , 


which  follows  from  the  second  part  of  (3.22)  . 


We  will  show  below  that  (3.22)  implies 


/ n • \ 

sequence  {q''  J } is  non-decreasing 


(3.25) 


ii'k 


Because  the  variance-covariance  matrix  Q must  have  a non-negative 

^ J 

determinant , 


r (2j)  2 . (2j ) (2j) 

Lq  . J < q^_  q2  * 


it  follows  that  the  sequence  in  (3.25)  is  ultimately  bounded  above  by  a 

g 0 

quantity  [ q^q2  3 2 + G where  e > 0 , and  so  this  sequence  must  have  a limit 


e r **  -s  c * * 1 

q . This  pertains  to  the  sequence  1^2 j > since  the  sequence  j_Q 2 j+i 

is  generated  by  the  same  type  of  recursion,  and  since  by  (3.24)  its  initial 


matrix  satisfies  the  analog  of  (3.22) , the  same  logic  proves  the 

existence  of  the  limiting  covariance  q° 

It  only  remains  to  prove  that  (3.22)  implies  (3.25).  Since  (3.22) 
implies  (3.23),  it  suffices  to  establish  the  prototype  induction  step, 
namely  that  (3.22)  implies 

> q • (3.26) 


Application  of  (3.8)  yields 

q(2)  = (a2+  S)“1[a4+  2Sa2+  (q^q)2]"1  {qa6 

+ [(q2”q)q-j+ A + q(s  + q1-q)]cr4+  [(q2-q)  A + q1(q2-q)2+  (s  + q1"q)A]a2 
+ (q2-q)2  A 3 

« (a2+  S)**1  [a4+  2Sa2  + (q2-q)2]“1  {qa6 
+ [2A  + qS]a4  + [2SA  + q^  (q2~q)2  ]a2+  (q^-q)^}  . 

Thus  (3.26)  is  equivalent  to 

qa6  + [2A  + qS  ]a4  + [2SA  + q^(q2~q)  ]a2  + (q2~q)2A 
> q(^2+  S)[a4  + 2Sa2  + (q2~q)2]  , 

or  in  turn  to 

(A  - qS){2a4  + [2S  + (q2~q)  ]a2  + (q2~q)2  ] > 0 , 


which  is  true  by  virtue  of  (3.18)  and  (3.22).  This  completes  the  proof 
that  the  limiting  covariances  (3.21)  exist. 
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matrix  satisfies  the  analog  of  (3.22) , the  same  logic  proves  the 

existence  of  the  limiting  covariance  q° 

It  only  remains  to  prove  that  (3.22)  implies  (3.25).  Since  (3.22) 
implies  (3.23),  it  suffices  to  establish  the  prototype  induction  step, 
namely  that  (3.22)  implies 

q(2)  > q . (3.26) 

Application  of  (3.8)  yields 

= (a2+  S)  1[c j4+  2Sa2+  (q2~q)2]  1 {q^6 

+ [(q2“ci)ci1+  A + <i(s  + q1~q)  t(q2-q)  a + q^(q2"q)2+  (s  + q1~q)^]cr 

+ (q2-q)2  a } 

= (a2+  S)"1  [a4+  2Sa2  + (q^q)2]"1  {qa6 
+ [2A  + qS]a4  + [2SA  + q1  (q2~q)2  ]a2+  (q2-q)2A}  » 

Thus  (3.26)  is  equivalent  to 

qa6  + [2A  + qS  ]a4  + [2SA  + q1(q2~q)2]a2  + (q2~q)2A 
> q(a2+  S)[ct4  + 2Sa2  + (q2-q)2]  , 

or  in  turn  to 

(A  - qS){2a4  + [2S  *f  (q2“q)  ]a2  + (q2~q)2  3 > 0 , 

which  is  true  by  virtue  of  (3.18)  and  (3.22).  This  completes  the  proof 
that  the  limiting  covariances  (3.21)  exist. 


- 

. 


Although  the  main  objectives  in  this  Section  have  been  achieved,  we 
should  note  several  items  of  "unfinished  business"  which  need  to  be 
settled  to  obtain  a thorough  analysis  : 

o e * 

(a)  Determine  whether  the  limiting  covariances  q and  q (and 

o e 

hence,  the  limiting  matrices  and  ) exist  in  general,  and  not  only 

under  the  extra  hypotheses  (3.22). 

(b)  Evaluate  the  various  limits  in  closed  form. 

(c)  If  (b)  is  too  ambitious,  at  least  determine  whether  the 
limits  are  independent  of  a 
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1.  Description  of  the  Algorithms  used  by  WHERSM 

The  algorithm  used  to  locate  each  unit  is  the  least  squares  squared 
(LSS)  algorithm.  Both  two-  and  three-dimensional  subroutines  are  included 
in  the  WHERSM  system.  The  three-dimensional  program  is  described  in 
detail  below. 

For  each  locator  i,  i = l,2,...,n,  let  (x^,  y , z ) be  its  estimated 
position,  r^  its  reported  distance  to  the  locatee,  and  w^  the  reciprocal 
of  its  variance,  the  calculation  of  which  is  described  below.  The  w.’s 
act  as  weights;  the  smaller  the  variance  of  a locator,  the  more  accurately 
we  know  its  position,  and  the  higher  its  weight  will  be.  Then  LSS  minimizes 


E = 


E.  = 

l 


n 


T 


i=l 


w.E. 

i l 


where 


7 


(1) 

(2) 


here  d_^  is  the  calculated  distance . between  (x^,  y , z ) and  (x,  y,  z)  , 
the  estimated  position  of  the  locatee,,  so 


d?  = (x-  xp2  + (y-  yi)2  + (z-  z^2 


(3) 


Note  that 


2 


2 


2 


2 


1 


I 
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When  the  estimated  positions  of  locators  and  locatee  are  near  their  true 

2 2 

positions,  r.  is  very  near  d.  , so  the  factor  (d.+  r.)  / 4r . is  very 
i i 2 1 1 1 

near  1 , and  behaves  very  much  like  (d^-r^)  * Thus  we  expect  this 

algorithm  to  produce  answers  nearly  identical  to  those  produced  by  the 

2 

least  squares  (LS)  algorithm,  for  which  E^  = (d^-  r^)  . Indeed,  computer 

tests  indicate  that  there  is  no  significant  difference  between  the  results 

of  these  two  algorithms.  Thus  LSS  is  used  because  it  is  more  economical, 

since  LS  requires  the  calculation  of  d^  , which  in  turn  requires  a 

square  root  to  be  taken  (see  (3)  ) . LSS  avoids  the  square  root,  since  only 

d?  is  needed, 
l 

Where  E is  a minimum,  its  partial  derivatives  with  respect  to  x,  y, 
and  z will  be  zero.  Thus  the  program  calculates 
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2 2 
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S ■ 2 (z 

i=l  r. 


Zi> 


(4) 


Since  d_^  depends  on  x,  y,  and  z (see  (3)  ),  F,  G,  and  H depend  non- 
linear ly  on  these  variables.  Then  the  system  (4)  cannot  be  solved  directly 
for  x,  y,  and  z . Instead,  LSS  does  successive  Newton-Raphson  iterations 
( [1],  pp.  201-223)  , beginning  with  a starting  point  .p  = (x,y,z)  . The 
first  iteration  yields  a new  point  p**  = (x  + Ax,  y +Ay,  z + Az)  ) this 
point  is  used  as  the  starting  point  for  the  next  iteration,  and  so  on.  The 
process  begins  with  the  first  order  approximation 


__  . bF  A . 3F  a . c>F  A 

F = F + r—  Ax  + — Ay  + — Az 
p p dx  dy  dz 


G * BG  . . 5G  a , dG  A 

p = G + — Ax  + — Ay  + — Az 

^ p Bx  dy  dz 


„ „ . dH  ' ,dH  A , dH  A 

H ^ = H + — Ax  +■  — Ay  + — Az 
p p dx  dy  dz 


Here  F denotes  the  value  of  F at  point  p , and  all  partial  derivatives 

are  evaluated  at  p . Thus  this  system  represents  the  expansions  of  the 

functions  F,  G,  and  H about  p . Now  we  want  F / , G '*  and  H > all  to 
9 P P P 

be  zero,  so  the  problem  reduces  to  that  of  finding  Ax,  Ay,  Az  which  satisfy 

the  three  simultaneous  equations 


|^Ax  + |^Ay  + — -Az+F  = 0 

dx  dy  dz  p 


— Ax  + Ay  + ~ Az  + G - 0 

dx  dy  dz  p 


(5) 


dH  A , dH  A , dH  A , ■ n 

— Ax  + T—  Ay+^-Az  + H = 0 . 

dx  dy  dz  p 


The  program  calculates  F,  G,  and  H (see  (4)  ),  and 
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then  solves  (5)  for  Ax,  Ay,  Az  . This  determines  p*=  (x  + Ax,  y + Ay, 
z + Az)  , which  is  used  as  the  starting  point  for  the  next  iteration.  For 
a given  set  of  locators,  the  program  will  do  up  to  4 iterations,  stopping 
if  an  iteration  converges,  namely 

n 

max(  | F | , | G | , | H | )<  S w , 

i=l 

or  if  an  iteration  diverges  (which  see  below)  . I.e.  convergence  occurs 
when  F,  G,  and  H are  all  small. 

Computer  tests  indicated  that  for  both  the  two-  and  three-dimensional 
versions  of  LSS  , 2 to  4 Newton-Raphson  iterations  are  generally  sufficient 
for  convergence,  even  when  the  starting  point  is  100  meters  from  the  true 
position. 

The  LSS  algorithm  allows  for  easy  elimination  of  outliers.  These 
are  locators  whose  reported  distance  (r^)  to  the  locatee  contains  a greater 
error  than  do  the  reported  distances  of  the  other  locators.  Let 

n n 

* ( r we)/  d w,  • 

i=i  11  i=i  1 


E 


(7) 


I 
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Now  any  locator  i is  eliminated  for  which 

E±>  4eT  . (8) 

Several  other  constants  greater  than  1 were  tried  in  place  of  4,  but 
computer  tests  showed  that  4 was  best. 

The  algorithm  iterates  on  all  locators  until  convergence  occurs , up 
to  4 iterations  (for  divergence,  see  below).  Then  E is  calculated  and 
some  locators  may  be  eliminated  via  (8)  . (If  no  locators  are  eliminated, 
we're  through).  Then  more  iterations  are  performed  on  the  reduced  set  of 
locators.  After  convergence,  (7)  and  (8)  are  again  applied  to  eliminate 
outliers  . This  iteration-elimination  process  continues  until  either  (1) 
no  locators  are  eliminated  in  a step^  (2)  1/4  of  the  locators  have  been 
eliminated,  or  (3)  the  algorithm  diverges  (see  below)  „ In  cases  (1)  and 
(2)  the  algorithm  terminates. 

There  are  two  types  of  divergence^  an  increase  in  the  objective  func- 
tion, E , from  one  iteration  to  the  next  and  a badly  scaled  determinant 
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less  than  10  in  absolute  value  which  makes  equations  (5)  impossible  to 
solve  accurately.  In  each  case,  we  halve  the  step  size  of  the  previous 
iteration,  and  try  again.  This  can  be  done  unless  the  iteration  at  which 
divergence  occurs  is  the  first  iteration  for  this  locatee . Then  we  must 
quit  without  finding  a solution. 

For  any  locatee  and  set  of  locators,  the  size  of  successive  iteration 
steps  should  decrease  rapidly.  If  not,  then  the  problem  is  ill-behaved, 
and  perhaps  dozens  of  iterations  will  be  necessary  to  achieve  convergence. 
To  avoid  this,  the  algorithm  attempts  to  " home  in"  on  the  solution 
faster.  If  the  step  size  in  the  x-direction  (or  y or  z)  is  Ax^  for  one 
iteration  and  Ax2  for  the  next,  and  if  | Ax^  | > % | Ax^ | , then  Ax2 

replaced  by 

2min(  |Ax^|,  | Ax2|  ) ® sgn  Ax2 


if  Ax^  and  Ax^  have  the  same  sign,  and  if  they  have  opposite  signs, 
AX2  is  replaced  by 


AX2  • Ax^  • 

Ax1 - Ax ^ 

A further  refinement  of  the  algorithm  was  tried;  this  involved  relocat- 
ing the  locators.  More  concerning  this  may  be  found  in  sec*  3 and  Working  Paper 
No.  8.  Once  the  Newton  iterations  have  converged,  there  is  an  error  between 
the  reported  distance  r_^  and  the  calculated  distance  d^  , for  each  locator 
i in  the  set.  By  changing  the  estimated  position  of  each  locator,  i.e.  by 
moving  it  along  the  straight  line  connecting  it  with  the  locatee,  this  error 
d - r_^  may  be  reduced,  even  eliminated.  After  relocating,  more  Newton- 
Raphson  iterations  may  be  done.  Several  runs  were  made,  using  different  values 
for  the  fraction,  f , of  the  error  eliminated  by  relocation,  with  0 < f < 1 . 
The  best  such  f was  the  smallest,  and  that  was  not  quite  as  good  as  no 
relocation  (f  = 0)  . Relocation  remains  an  option  in  the  LSS  routines,  but 
one  which  is  not  now  used. 

The  2-dimensional  LSS  program  is  analogous  to  the  3-dimensional  one.  It 
solves  the  equations 


dF 

dx 


Ax  + 


dF 

dy 


Ay  + F **  0 


dG 

dx 


Ax  + ~ Ay  + 

dy 


G = 0 


Here,  F and  G and  their  partial  derivatives  are  as  described  above,  except 
that  all  z-components  are  zero.  The  2-dimensional  program  is  less  likely  to 
result  in  divergence  of  an  iteration,  so  the  special  steps  to  handle  divergence 
(see  previous  page)  are  omitted. 


The  weights  w (see  (1))  are  again  taken  as  the  reciprocals  of  the 
2 1 

variances  o\  . When  a unit  has  been  located,  the  variance  of  its  location 
is  calculated  as 


n 


( £ w.E.) 

- l l 
i=l 


n 


(n-m)  £ w . ) 

i=l  1 


C9) 


when  n is  the  number  of  locators,  m the  dimensionality  (2  or  3)  of  the 
algorithm* 


When  the  locatee  is  lost,  so  that  no  starting  point  exists  for  the  LSS 
routine,  the  least  squares  linear  (LSL)  routine  is  called  to  provide  a start- 
ing point.  First,  a slant-range  reduction  is  done  to  convert  3-D  reported 
distances  to  2-B  . Then  the  2-D  LSL  algorithm  is  called.  This  is  the  same 
basic  algorithm  which  is  described  in  sec.  2 of  [2]  * The  modifications  are 
these:  all  pairs  of  circles  are  used,  and  outliers  are  eliminated.  Consider 
the  M radical  axis  lines  each  formed  by  a pair  of  circles ) the  equation 
for  such  a line  is  found  by  subtracting  the  equation  of  one  circle  from  the 
equation  of  the  other  as  in  [2]  , p » 10„  For  each  such  line  k,  k — 1,2, ...M, 

let  d,  be  the  calculated  distance  from  the  locatee  to  the  line.  Then  LSL 
k 

minimizes 


D 


M 

L 

k=l 


After  a locatee  has  been  found  using  all  M lines,  the  routine  calculates 


d2  = D/M  , 

the  average  squared  distance  to  the  lines.  Then  any  line  k is  eliminated 
for  which 


j- 


The  reduced  set  of  locators  is  used  to  calculate  a position,  d is 
recalculated,  and  a line  k is  eliminated  for  which 


2 2 
df>  4 • dZ 
k 


The  process  is  repeated  twice  more,  once  with  the  elimination  condition 


d.  > 2 

k 


2 2 
df  > d 
k 


and  finally  with 


The  LSL  algorithm  will  fail  only  when  the  equations  at  the  bottom  of 
p.  13  of  [2]  are  ill-conditioned.  This  will  happen  when  the  M lines  are 
nearly  parallel,  i0e0  the  locators  all  lie  near  one  straight  line.  In  this 
case,  no  starting  point  can  be  found  for  LSS 


The  routine  FINDIT  calls  LSS  (2-D  or  3-D)  and  LSL  when  necessary. 

If  the  locatee  is  not  lost,  two  previous  locations  are  used  to  extrapolate 
in  the  x-  and  y-  directions  to  get  a starting  point  for  LSS.  The  last 
z-value  is  used,  with  no  extrapolation.  Then  LSS  is  called;  the  2-dimen- 
sional  routine  for  ground  units,  3-dimensional  for  airplanes. 

If  the  unit  is  lost,  LSL  is  called  after  doing  a slant-range  reduction 
(u  sing  the  last  known  z-value)  to  convert  3-dimensional  distances  to  planar 
distances.  Then  using  the  planar  distances,  the  2-dimensional  LSS  is 
called.  If  the  unit  is  an  airplane,  the  x-  and  y-coordinates  returned  by 
LSS  (2-D)  plus  the  last  known  z-value  are  used  as  a starting  point  for 
LSS (3-D) . 

Control  is  returned  to  FINDIT  with  the  variance  defined  by  (9)  . This 
represents  the  variance  of  the  position  estimate,  based  on  the  reported 


. 


distances.  The  total  variance  is  computed  by  adding  RVAR,  the  variance 
of  the  distances.  If  the  locatee  is  known  to  be  moving,  its  weight  is 
calculated  by  inverting  the  total  variance. 


If  the  unit  is  fixed,  the  current  position  estimate  is  taken  to  be 
a weighted  average  of  the  old  estimate  and  the  just-calculated  estimate. 
Then  the  weight  of  the  locatee  is  computed  as  a weighted  average  of  the 
old  weight  and  the  just-calculated  weight.  A further  description  of 
these  calculations  may  be  found  in  Section  three  of  this  report. 

2.  Description  of  Algorithms  Programmed  for  WHERSM  But  Not  Now  Used 

The  linear  method  (LM)  and  smallest  tangent  circle  (STC)  method 
are  described  in  [2],  They  were  both  found  to  be  less  accurate  than 
the  LSS  and  LSL  algorithms. 

The  least  squares  method  mentioned  in  §1  is  the  same  as  LSS, 
except  equations  (2),  (4),  and  (6)  are  replaced  by 
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/ 
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The  other  algorithm  programmed  was  MINMAX  which  minimizes  the  maximum 
of  the  errors  d^  - r'  . This  algorithm  proved 

to  be  less  accurate  than  STC  and  much  less  accurate  than  LSS.  Further, 
it  uses  10  to  100  times  the  computer  time  the  other  algorithms  use, 
since  (in  the  two-dimensional  case)  it  considers  all  pairs  and  triples 


of  locator  circles. 
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ABSTRACT 


In  Working  Paper  No.  2 it  was  suggested  that  since  the  locators'  posi- 
tions are  not  known  exactly,  it  might  prove  profitable  to  re-estimate  these 
positions  as  an  adjunct  to  estimating  the  position  of  the  "locatee"  unit. 
This  paper  presents  an  iterative  method,  along  the  lines  proposed  in  Working 
Paper  No.  2,  for  handling  the  computations  of  such  a "relocation"  process. 
The  discussion  is  carried  out  in  the  context  of  the  LS(least  squares)  and 
LSS(least  squares,  squared)  position-location  methods.  As  background, 
evidence  is  given  which  shows  that  these  two  criteria  behave  very  similarly 
so  long  as  range  measurement  errors  remain  small. 


Note : Project  working  papers  are  informal  documents  prepared  to  facilitate 

discussion  and  communication;  they  may  contain  tentative  or  relatively 
unchecked  material. 


CONTENTS 


0o  INTRODUCTION - 1 

1.  RELATION  BETWEEN  LS  AND  LSS  APPROACHES 3 

2.  THE  RELOCATION  PROCESS - 7 


0.  'INTRODUCTION 


In  this  paper  the  term  location-operation  will  be  used  for  a process 
of  attributing  a position  to  one  unit  (the  locatee)  on  the  basis  of  its 
recorded  distances  from  n other  units  (the  locators)  . These  recorded 
distances  are  of  course  subject  to  measurement  error. 

Roughly  speaking , the  situation  studied  in  Project  "WHERE"  is  one  in 
which  the  role  of  "locatee"  varies  over  some  totality  of  units , with  many 
of  these  units  in  motion  at  any  time.  From  this  it  is  apparent  that  for 
any  one  location-operation,  the  positions  of  the  locators  (which  may  have 
been  the  locatees  of  previous  operations) are  not  in  general  known  exactly. 
This  point  was  emphasized  in  Section  5 of  Working  Paper  No.  2 ["  An  Itera- 
tive Approach  to  the  Location-Operation"  , APJ.  Goldman,  10/71],  where  it 
was  suggested  that  the  location-operation  be  expanded  to  include  re-estima- 
tion  (i.e.,  relocation)  of  the  locators'  positions,  based  on  the  latest  set 
of  range  measurements. 

The  present  paper  has  as  its  main  purpose  the  computational  specifica- 
tion of  such  a relocation  procedure.  This  is  given  in  Section  2,  which  also 
discusses  the  respective  merits  of  the  LS  (least  squares)  and  LSS  (least 
squares,  squared)  position-location  methods  in  this  context.  Section  1 pre- 
sents some  relevant  background  information  concerning  the  relationship 
between  these  two  algorithms. 

We  conclude  this  Introduction  by  introducing  some  necessary  mathematical 
notation: 

0 •»  (arbitrary)  origin  of  coordinate  system  , 

X **  position  to  be  atrributed  to  locatee, 
x = vector  OX 

A^  *=  position  to  be  attributed  to  i-th  observer 
(i  1,  2,o.e,  n) 


> ‘ 


* 


' 


vector  OA. 

1 


a. 

1 


initial  estimate  of  a. 

1 


— |jx  - a^  II  55  estimated  distance  between  locatee 

and  i-th  observer , 


P 


i 


measured  value  of  this  distance  . 


If  the  distinction  between  a^  and  is  dropped  (i.e.,  the  initial 

estimates  of  the  locators'  position  are  accepted  without  revision),  then  we 
have  the  "pure"  location-operation,  which  can  be  expressed  as  the  problem  of 
making  a "good"  choice  of  x , given  the  and  a^  . The  LS  method 

expresses  this  problem  as  that  of  minimizing  the  quadratic  "penalty  function" 

f(x)  = D w^(rj,  - pp2  (0.1) 

i 

where  the  w^  are  positive  "weights"  and  x enters  via  the  r^s  $ 
the  LSS  approach  involves  the  minimization  of 

g(x)  = S - P?)2  , (0,2) 

with  "weights"  w^  . The  modification  of  these  formulations,  to  encompass 
relocation,  will  be  taken  up  in  Section  2 . 


1„  RELATION  BETWEEN  LS  AND  LSS  APPROACHES 


In  this  section,  as  background  for  the  discussion  of  relocation  in 
Section  2,  we  discuss  the  relation  between  the  LS  and  LSS  approaches 
to  the  "pure"  location-operation.  Three  pieces  of  evidence  will  be  given 

to  show  that  if  the  "weights"  in  these  two  criteria  the  w^  of  (0*1) 

and  the  w^  of  (0o2)  are  properly  related,  then  the  two  methods 

will  not  differ  significantly  in  accuracy  so  long  as  the  tracking  process 
is  under  control 0 (The  last  phrase  means  that  the  units'  positions  remain 
well  enough  known  that  the  discrepancies  between  measured  ranges  p_^  and 
calculated  ranges  r^  are  small  relative  to  the  magnitudes  of  the  ranges 
themselves .) 

The  first  piece  of  evidence  begins  with  the  observation  that  the 
"under  control"  assumption,  r^/p^  ~ 1 > yields 


,2  2.2  , , .2  , .2 
<ri  - p.)  ■=  <r  + PjL)  (r.  - p.) 


(1  + r^p.)2  PJ(r.  - p.)2 


2 2 
4p7(r . -r  p.)  , 

is  l rx  9 


or  equivalently 


2 2 2 2 
(l/4p2)(r^  - ppZ 


(ri  " pi> 


2 * 

Thus,  if  we  set  w^  « (l/4p^)  w^  , then  comparison  of  (0 0 1 ) with  (0o2)  shows 
that  g(x)  ex  f(x)  ; equivalently,  if  we  set 


w 


i 


2 

i 


(1.1) 


then  comparison  of  ( 0 0 1 ) yields-  g(x)  « 4f(x)  , so  that  minimization  of 
f(x)---  i.e.,  the  LS  approach  ---  is  "approximately  equivalent"  to  the 

LSS  approach  (minimization  of  g(x)  ) . 

Equation  (1*1)  gives  the  appropriate  relation  between  the  two  sets  of 
weights.  As  yet,  no  rationale  has  been  found  for  assigning  to  the  w*  any 
value  other  than  unity.  In  fact,  since  the  distribution  of  errors  in 
(unrejected)  ranges  seems  (by  best  engineering  judgement  at  this  time)  to 
be  Gaussian,  with  range- independent  variance  and  zero  mean  (except  for  a 
small  percentage  of  ranges  biased  by  about  8m),  maximum  likelihood  estima- 
tion would  strongly  recommend^^  use  of  LS  with  w^  - 1 (together  with 
some  technique  for  removing  "outliers"  ) , which  by  the  above  is  approximately 
equivalent  to  use  of  LSS  with 

w.  = 1/p*  . (1.2) 

This  choice  of  w*  and  w^  will  be  used  throughout  the  rest  of  this  paper. 

The  preceding  analysis  shows  that  the  penalty  functions  for  LS  and 
LSS  are  approximately  equal  (when  measured  and  calculated  ranges  are  close 
together),  but  one  would  like  some  direct  assurance  that  the  position-estimates 
obtained  by  minimizing  these  two  functions  are  also  nearly  equal.  This  is 
provided  by  our  second  piece  of  evidence,  an  analysis  for  a simple  1-dimen- 
sional  configuration  in  which  the  locators  are  situated  at  the  points  with 
coordinates  R and  (~R)  with  R large , while  the  locatee  is  actually 
situated  at  the  origin.  The  LS  approach  calls  for  the  minimization  of 

f(x)  = [(R-x)  - Q1  ]2  + [(R-fx)  - P2]2 

2 

= 2 {x  + (p^-  p^)x  } + const.. 

See  Section  1 of  Working  Paper  No.  2 * 


(1) 


. 


which  is  readily  found  to  be  given  by 


x = (P2»  p ^ ) /2  . (LS  solution)  (103) 

The  LSS  approach  calls  for  the  minimization  of 

/ x -2r/T.  .2  2.2,  -2r/„ , .2  2.2 

g(x)  = Px  [(R-x)  - Px  ] + p2  [(R+x)  - p2] 

-2  4 -2  4 2 

- p^  (R-x)  + P2  (R+x)  - 4x  + const,  j 

setting  dg/dx  = 0 leads  to  the  cubic  equation 

, 2 , 2.  3 OT>/  2 2.  2 ro  2,  2,2.  0 2 2 . , 3,  2 2 . „ 

(P1  + P2)x  + 3R(pL“  P2)x  + [3R  (Px  + P2)  - 2pxP2  ] x + R (p1»  P2  ) * 0 . (1 

Let  E denote  a measure  of  maximum  range-error.  In  the  worst-case  scenario 
described  by  p^  = R-E  and  p2  = R + E , the  LS  and  LSS  approaches  yield 

exactly  the  same  answer,  both  (1„3)  and  (1„4)  giving  x « E . In  the  alter- 
native scenario  described  by  p^  = R and  p2  = R + E , the  respective  solutions 
are 

*LS  = !E  » "LSS  - IE  [1  - |(E/R)2  ] (1.5) 

where  the  expression  for  x^gg  omits  higher-order  terms  in  the  small  quantity 
E/R  ; if  for  example  R is  at  least' 1 km  and  E is  at  most  20m  , then  the 
two  position-estimates  differ  by  at  most  0o0015  m . 


Our  third  piece  of  evidence  is  empirical  in  nature.  For  each  of  two 
distributions  of  range -measurement  errors  (’’small"  and  "large",  respectively) 
72  trial  2-dimensional  location-operations  were  carried  out  using  both  LS 
and  LSS  . Each  of  18  units  served  as  locatee  4 times,  (for  each  error  dis- 
tribution and  each  method),  with  the  other  17  units  serving  as  locators.  For 
the  "small"  case,  the  locations  were  as  much  as  3.32m  in  error,  but  the  LS 
and  LSS  results  differed  from  each  other  by  no  more  than  0„ 0004m  ; for  the 
"large"  case,  despite  location  errors  up  to  31.9m,  the  two  methods  yielded 
locations  no  more  than  0o045m  apart. 


2.  THE  RELOCATION  PROCESS 


When  the  location-operation  is  enlarged  to  include  relocation  (of  the 
locators) , it  can  be  formulated  as  making  a ’‘good"  choice  of  x and  the 
a^  , given  the  and  . The  approach  given  in  Section  4 of  Working 

Paper  No.  2 is  based  on  two  notions,  one  of  which  we  have  found  it  is  useful 
to  modify  in  a manner  described  below. 

The  first  notion  is  that  the  choices  of  x and  a!s  should  be  made  so 

i 

as  to  minimize  a penalty  function  which  generalizes  those  given  earlier  [see 
(0.1)  and  (0.2)]  for  "pure"  location.  Specifically,  this  function  is 


f (x$  a1,...,  an)  = 


2.  w'(r.  - p.)^  + D.  W.  I!  a.-  o'.  I!  (4.1a) 

iiNi  l i l 11  l i ■'  v ' 


for  the  LS  approach,  and 


g(x;  , . . . , a^) 


2 2 2 

£.  w.(rT  - PV  ) + £.  W. 

l iNi  i li 


(4.1b) 


for  the  LSS  approach  ; the  are  positive  numerical  weights. 

The  second  notion  is  to  accomplish  this  minimization  by  a process 
which  alternates  between  location  steps  (in  which  the  a^  are  held  fixed 
and  an  appropriate  x is  chosen)  and  relocation  steps  (in  which  x is  held 
fixed  and  appropriate  a^'s  are  chosen.  This  process  can  be  formalized  as 
follows : 

INITIALIZATION:  Set  k « 1 , and  = Oi . (all  i ) . 

9 l l ' 

(k) 

STEP  1 (location,  k-th  pass)  : Given  the  current  estimates  a; 

(k)  1 

of  the  locator  positions,  choose  x - xv  ' to  minimize  the  first  sum  in 

(4.1),  i.e.  the  part  of  (4.1)  which  depends  on  x 


STEP  2 (relocation,  k-th  pass)  : For  each  i,  keep  x fixed  at 

(\r\ 

xv  ' and  choose  a.  to  minimize  the  sum  of  those  terms  in  (4.1)  which 

1 (k) 

depend  on  a^  . Increase  k by  1 , and  set  each  a^  J equal  to  the 

a_^  just  found.  Return  to  Step  1 . 

For  Step  1,  the  material  in  Section  1 shows  that  there  is  little 
reason  to  prefer  one  of  LS  or  LSS  over  the  other  so  far  as  location- 
errors  are  concerned.  We  therefore  select  LSS,  since  it  seems  likely 
to  converge  more  quickly  and  requires  fewer  time-consuming  calculations 
per  iteration.  Thus  Step  1 employs  (4.1b)  rather  than  (4.1a),  and  so 
according  to  the  description  above,  the  treatment  of  the  i-th  locator  in 
the  k-th  relocation  step  would  require  choosing  a^  to  minimize 


w. 

l 


- a . 

l 


(4.2) 


For  Step  2,  however,  it  is  more  convenient  to  minimize  instead  the  sum  of 
corresponding  summands  from  (4.1a_),  namely 


Hap 


(4.3) 


The  over-all  method  becomes  a kind  of  LS-LSS  hybrid  ; Steps  1 and  2 in 
effect  involve  different  penalty  functions,  so  that  the  "first  notion" 
described  above  has  indeed  been  modified. 

In  Working  Paper  No.  2 it  was  envisiged  that  the  alternating  process 
would  be  iterated  until  it  had  "settled  down"  satisfactorily.  This  idea, 
too,  has  been  modified  after  further  reflection.  The  relocation  of  the 
i-th  locator  is  based  mainly  upon  a single  new  datum,  namely  ; it  seems 


unwise  to  put  too  much  reliance  -or  to  base  too  much  calculation  on  this  one 
datum  (which  is  subject  to  measurement-error).  Besides,  from  the  viewpoint 
of  reducing  total  computation  time  it  is  desirable  to  limit  the  number  of 
iterations.  Thus  a sensible  compromise  seems  to  be:  locate  the  locatee 

(Step  1);  relocate  the  locators  (Step  2)  ; then  locate  the  locatee  again 
(second  pass  through  Step  1),  and  stop.  This  truncated  version  of  the  full 
process  is  assumed  below. 

We  turn  now  to  the  execution  of  Step  2.  With  superscripts  (k)  and 
subscripts  (i)  dropped  for  simplicity,  the  function  (4.3)  to  be  minimized 
can  be  written 

h(a)  = w'{|jx-a  ||  - p } ^ + W ||  a-a  ||2  . (4.4) 


Consider  any  trial  position  a'  of  a .It  lies  on  a sphere  of  radius 
||  x-a^  ||  centered  at  x . Suppose  a * is  replaced  by  the  point  of  that 
same  sphere  which  is  closest  to  (X  ; then  the  first  summand  in  (4.4)  is 
unchanged,  while  the  second  summand  is  reduced-  Hence  the  minimum  can  only 
occur  for  such  a "closest  point".  Simple  geometric  considerations  show 
that  such  a point  lies  either  on  the  line  segment  [x, <x  ] or  on  its  extension 
through  a . Thus,  for  some  scalar  s > 0 , we  have 

a = x + s(o?  - x)  , 

leading  to 

x - a = s ( x-a  ) , a - cx  = (1-s)  (x-o)  . 


Substitution  of  these  results  into  (4.4)  converts  h(a)  into  a function 
of  s \ with  the  abbreviation 


r 


I!  x - a 


(4.5) 


this  function  reads 


H(s)  = w'{sr-p)^  + Wr^(l-s)^ 


= (w^d-  W)r^s^  - 2r(w^p  + Wr)s  + (w^p^  + Wr^  ) , 


(4.6) 


and  is  minimized  by  setting 


s *»  (w^p  + Wr)/(w**  + W)r  . (4.7) 

(The  availability  of  such  a simple  explicit  solution  to  the  minimization 
problem  justifies  our  earlier  assertion  that  use  of  (4.4)  rather  than  (4.2) 
is  more  convenient.)  If  p > r then  s > 1 , signifying  that  a lies  on 
the  extension  of  segment  [x,  O']  through  a (i.e.,  relocation  has  moved 
the  locator  farther  from  the  locatee)  ; if  p < r , then  s < 1 and  a 
lies  on  the  segment  (so  that  the  locator  has  been  moved  closer  to  the  locatee). 
These  results,  which  do  not  depend  on  w and  W ',  seem  quite  consistent  with 
what  one  would  expect  intuitively. 

The  only  remaining  question  is  the  selection  of  numerical  values  for  the 
weights  W^  . We  have  at  present  no  definitive  answer  to  this  problem. 

There  are  two  distinct  considerations  to  be  taken  into  account: 

(a)  If  we  consider  the  current  location  (and  relocation)  operation  in 
isolation,  then  the  weight  W_^  should  depend  primarily  on  our  confidence  in 
the  accuracy  of  the  initial  estimate  Oi ^ of  the  i-th  locator's  position 
(relative  to  the  accuracies  of  the  other  ct  ' s , and  also  of  the  range  measure- 
ments). In  an  idealized  Gaussian-error  case,  for  example,  max imum- likelihood 
estimation  would  dictate  weights  inversely  proportional  to  the  respective 
variances  of  the  associated  random  errors.  For  our  practical  situation,  how 
can  we  quantify  the  degree  of  confidence  merited  by  the  initial  estimate  Q\^  ? 


It  may  well  prove  possible  to  concoct  an  appropriate  measure  based  on  the 
sizes  of  the  discrepancies  | r_.  - p | obtained  during  the  last  previous 
location-operation  in  which  the  i-th  of  the  current  locators  figured  as 

locatee probably  modified  by  some  "decay  factor"  to  take  account  of  how 

much  time  has  elapsed  since  this  earlier  operation. 


(b)  However , one  needs  to  consider  the  "current"  operation  in  the 
light  of  the  need  to  maintain  stability  of  the  entire  location  process;  the 
selection  (or  least  the  effect)  of  the  must  reflect  not  only  the  one- 

shot  independent -measurement  case  as  in  (a),  but  also  the  sequential  corre- 
lated-error  many-locator  nature  of  the  real  problem.  If  a relocat ion(based 
mainly  on  a single  and  possible  "bad"  range  measurement  ) is  allowed 

to  be  too  large,  its  influence  would  persist  through  a number  of  operations 
until  the  unit  in  question  again  took  its  turn  as  locatee,  and  so  such 
instabilities  in  estimation  of  positions  could  result  in  large  errors. 

The  following  tactic  may  prove  an  acceptable  compromise.  Weights 

will  be  selected  according  to  the  considerations  in  (a)  above.  The 

resultant  solutions  (4.7)  of  the  n minimization  problems  for  relocation 

(1  problem  for  each  of  the  n locators)  will  dictate  moving  the  estimated 

position  of  the  i-th  locator  a certain  distance  (from  a.  ) in  a certain 

(k)  1 

direction  along  the  line  through  xv  7 and  & . To  "temper"  this  solu- 
tion in  accordance  with  the  considerations  in  (b),  we  will  "damp"  that 
distance  by  a multiplicative  factor  (say,  between  0.1  and  0.6)  before 
actually  revising  the  estimate  of  the  locator's  position.  Numerical  exper- 
imentation will  be  required  to  test  the  performance  of  this  approach  and  to 
determine  what  values  of  the  adjustment  factor  are  preferable. 
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1.  INTRODUCTION 


(*) 

WHERSM  is  a computer  code  that  simulates  the  operation  of  a 
Micro-Navigation  Position  Location  System  (MNPLS) . It  provides  a "real 
life"  framework  within  which  the  propagation  over  time  of  position  loca- 
tion errors,  resulting  from  the  operation  of  such  a system,  can  be  deter- 
mined. It  further  provides  a facility  for  monitoring  the  operation  of 
such  a system  under  a variety  of  movement  scenarios,  and  using  a variety 
of  algorithms  for  accomplishing  the  location  operation. 

This  paper  documents  the  present  version  of  the  WHERSM  code  at  three 
levels:  the  "executive"  level,  the  "user's"  level,  and  the  "analyst's" 

level.  This  introductory  section,  which  serves  as  the  executive  level 
of  documentation,  is  aimed  at  the  person  who  is  interested  only  in  an 
overview  of  what  WHERSM  is  and  what  it  does.  The  user's  level  of  docu- 
mentation is  meant  for  one  who  is  concerned  with  making  a computer  run 
with  the  code  as  it  now  stands;  i.e.,  with  no  major  program  modifications. 
Sections  10  and  11  provide  documentation  of  all  the  inputs  and  outputs  — 
the  information  needed  to  make  a run.  The  remainder  of  the  paper, 

Sections  2 through  9 and  appendices  A,  B,  and  C,  is  designed  with  the 
analyst  in  mind.  By  an  "analyst"  is  meant  one  who  needs  not  only  to  run 
the  program,  but  also  perhaps  to  modify  it  to  handle  somewhat  different 
situations.  This  paper  should  include  enough  information  to  start  him 
on  his  way. 

Due  to  the  nature  of  the  MNPLS  under  study,  WHERSM  is  a "fixed  time 
increment"  simulation;  i.e.,  simulated  time  is  advanced  by  a constant 
increment.  At,  rather  than  by  variable  increments  dictated  by  the 
occurrence  of  events,  as  would  be  the  case  in  the  type  of  simulation 
classified  as  "next  event". 


(*) 


Full  name:  WHERESIM.  Compressed  because  of  computer  restrictions 

to  6 characters  per  name. 


The  program  is  written  in  FORTRAN  V,  but  with  very  little  conversion 
it  can  be  compiled  with  a FORTRAN  IV  compiler.  It  is  liberally  annotated 
with  comments,  and  although  the  code  documented  here  employs  program  over- 
lays, parallel  processing,  and  drum  usage,  every  attempt  has  been  made  to 
keep  machine  dependence  minimized.  Furthermore,  the  WHERSM  code  is  modular 
in  nature,  so  that  the  logic  flow  is  kept  simple  and  easy  to  understand, 
and  so  that  different  concepts  in  deployment,  movement,  intervisibility, 
position  location  method,  etc.,  can  with  a minimum  of  reprogramming  be 
implemented  and  tested.  This  modularity  also  allows  extreme  flexibility 
in  the  type  and  quantity  of  inputs  required  to  make  a run. 

The  basic  ideas  involved  in  an  implementation  of  the  MNPLS  are  that 
there  are  n units  in  the  field,  any  of  which  at  each  moment  might  (in 
the  most  general  case)  be  either  moving  or  fixed.  Every  At  seconds,  a 
master  computer  receives  information  in  the  form  of  ranges,  subject  to 
measurement  error,  from  one  of  these  units  to  some  or  all  of  the  others. 

The  master  computer  must  then,  using  position  location  algorithms,  estimate 
the  position  of  that  unit  and  also  prepare  to  receive  the  set  of  ranges 
for  the  next  unit  in  sequence.  The  elements  of  this  description  that 
are  relevant  to  a simulation  study  of  an  MNPLS  are:  movement,  range 

computation  (including  the  occurrence  of  measurement  error)  , and  position 
estimation. 

In  particular,  the.  simulation  must  have  some  facility  for  determining 
the  true  positions  of  all  the  units  at  all  relevant  times.  It  Is  easy  to 
provide  a movement  scenario  in  terms  of  movements  in  the  xy-plane  (longi- 
tude-latitude) . A starting  x and  y,  an  azimuth  of  movement,  and  a 

(*) 

speed  will  suffice  for  this  . It  is  the  z coordinate,  the  height  above 
sea  level,  that  causes  difficulty. 


(*)  As  will  be  explained  in  Section  2,  the  movements  are  assumed  to  be 
piecewise  linear,  so  that  in  fact  a set  of  azimuths,  speeds,  and 
change  times  are  required. 


■ 


■ 


There  are  at  least  two  methods  of  providing  the  z coordinate  in 

a simulation.  One  involves  a continuous  function  which  yields  a value 

of  z for  any  values  of  x and  y;  the  other  is  by  means  of  a discrete 

function  which  specifies  values  of  z for  particular  values  of  x and 

y.  In  the  discrete  case,  values  of  x and  y other  than  the  prescribed 

ones  are  assumed  to  coincide  with  the  nearest  available  grid  point,  whose 

(*) 

specified  height  is  then  used  . Both  of  these  techniques  are  explained 
in  detail  in  section  2. 

The  next  item  with  which  the  simulation  of  the  MNPLS  must  be  concerned 
is  the  computation  of  ranges.  Of  course,  given  the  ability  to  compute  the 
true  positions  of  all  units  at  any  instant  of  time,  the  true  ranges  can 
easily  be  calculated,  since  they  are  simply  Euclidean  distances.  But,  in 
the  "real  life"  situation,  true  ranges  are  never  reported  to  the  master 
computer  due  to  measurement  errors  intrinsic  to  the  particular  method  of 
measuring  those  ranges.  This  phenomenon  is  easily  simulated  once  a 
probability  distribution  of  measurement  errors  is  specified,  since  the 
simulation  program  can  simply  compute  pseudo-random  numbers  which  satisfy 
such  a distribution  law. 

Another  point  relevant  to  range  computation  is  that  a simulation  must 

he  capable  of  determining  whether  a range  should  be  computed;  i.e., 

whether,  for  a specified  pair  of  units  at  a specific  moment,  a radio  line 

(**) 

of  sight  (LOS)  exists.  There  are  many  ways  of  providing  a simulation 

program  with  the  ability  to  determine  the  existence  or  non-existence  of 


these  intervisibilities.  The  method  used  in  WHERSM  is  to  maintain  a set 


(***) 


of  "intervisibility  snapshot  matrices",  I..  (i  = l,2,...,n;  j = l,2,...,n; 

i J t • 


(*)  Four-point  interpolation  could  have  been  used  here,  but  is  very  time 
consuming.  And  since  it  was  determined  that  no  significant  further 
loss  of  accuracy  occurs  with  the  "nearest  point"  method,  that  one 
was  implemented. 

(**)  See  p.28  of  "First  Interim  Progress  Report  on  NBS  Project  WHERE  in 
Support  of  USAAMCA  Project  MNPLS",  NBS  Report  10663,’  December  1971. 

(***)  Obtained  from  the  Electromagnetic  Compatibility  Analysis  Center, 
Annapolis,  Md. 


■ 
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t = 0,.  20,  40,  ...,  120)  where'  I.  = 1 if  unit  j can  receive  and  report 

•t  J t 

a ranging  signal  from  unit  i at  time  t minutes  into  the  simulation,  and 

I.  = 0 if  it  cannot.  Furthermore,  in  order  to  avoid  drastic  changes 
ijt 

in  the  intervisibility  pattern  at  the  discrete  times  indicated  by  values 
of  the  subscript  t,  each  row  of  a particular  intervisibility  matrix  defined 
by  a specific  value  of  t,  say  is  overlaid  with  the  corresponding  row 

of  the  next-in-sequence  matrix  (defined  by  t'  + 20)  at  some  randomly  chosen 
time  between  t'  and  t'  +20.  (These  times  are  chosen  independently 
for  each  row  and  for  each  value  of  t.)  More  detail  regarding  each  aspect 
of  range  computation  will  be  found  in  sections  3 and  5. 

The  last  item  of  concern  for  the  MNPLS  simulation  is  that  of  position 
estimation.  Although  that  item  is  an  integral  part  of  any  simulation  study 
of  this  nature  (indeed,  the  most  important  part),  no  discussion  of  the 
various  algorithms  used  in  WHERSM  is  provided  here,  since  the  project  staff 
members  responsible  for  the  development  and  implementation  of  position 
location  algorithms  have  presented  such  a discussion  in  Working  Paper 
Number  7. 

In  order  to  provide  the  reader  with  a better  understanding  of  the 
flow  of  events  within  WHERSM,  and  to  facilitate  further  discussion  of  the 
WHERSM  program,  a flowchart  is  given  in  figure  1 which  indicates  the  basic 
logic  of  the  simulation.  With  one  exception  (see  section  9),  each  box  in 
the  flowchart  corresponds  to  an  easily  recognizable  module  within  the 
WHERSM  code  listed  in  appendix  B.  In  every  case,  entry  to  box  2 is  made 
just  after  the  simulation’s  current  time  is  incremented  by  a constant 
value.  That  constant.  At,  is  referred  to  as  the  sub cycle  time  and  is  the 
time  during  which  all  the  operations  required  to  locate  a given  unit  Cone 
completion  of  the  loop  formed  by  boxes  2 through  7 of  the  flowchart)  must 
be  performed.  Cycle  time  is  the  time  required  to  locate  each  unit  a+  least 
once.  This  may  be  larger  than  n • At,  where  n is  the  total  number  of 
units,  since  certain  faster  moving  units  may  need  to  be  located  more  often 


than  do  slow  moving  units.  Locatee  refers  to  that  unit  which,  at  any 
instant  of  simulated  time,  is  in  the  process  of  being  located,  and 

locator  refers  to  a unit  whose  range  to  the  locatee  is  used  to  estimate 
the  position  of  the  locatee. 

The  remainder  of  this  report  discusses  the  workings of  WHERSM  at  a 
more  detailed  level.  Section  2 explains  Boxes  2 and  3 of  the  flowchart 
m figure  1,  while  Section  3 discusses  Box  4.  Section  4 explains  Boxes 
5,  6,  and  7;  Section  5,  Boxes  8 and  9;  Section  6,  Boxes  10  and  11; 

Section  7,  Box  12,  Section  8,  Boxes  13  and  14;  and  Section  9,  Box  15. 
Section  10  contains  a discussion  of  the  input  needed  to  make  a run  with 
the  present  version  of  the  code;  and  Section  11  contains  a discussion, 
with  examples,  of  the  outputs  that  will  result  from  such  a run.  Appendix 
A is  a dictionary  of  variables  used  in  the  program;  Appendix  B,  a listing 
of  the  program;  and  Appendix  C,  a very  detailed  flowchart  of  the  program. 

Although  the  remaining  section  headings  refer  directly  to  the  flow- 
chart in  figure  1,  the  content  of  those  sections  can  be  better  understood 
by  referring  to  Appendices  B and  C. 


. 


2.  BOXES  2 and  3;  THE  MOVEMENT  MODULES 


The  movement  modules  are  responsible  for  computing  the  true  positions 
of  each  unit  every  At  seconds.  In  the  current  version  of  WHERSM,  this 
is  a two-stage  process;  i.e.,  x and  y are  computed,  and  then  z is 
computed  using  x and  y.  The  method  by  which  x and  y are  computed 
is  explained  first. 

In  the  xy-plane,  units  are  assumed  to  move  in  a piecewise  linear 
manner  (the  projections  of  the  three-dimensional  curves  of  motion  onto 
the  xy-plane  are  piecewise  linear) . The  piecewise  linear  motion  is 
effected  by  maintaining  a matrix  of  current  true  positions  (TP)  and  a 
matrix  of  current  velocity  vectors  (VV) , both  of  which  are  initialized 
by  card  input,  and  by  updating  TP  linearly  as  follows: 

TP (1,1)  = TP  (1,1)  + W (1,2) 

TP  (1,2)  = TP  (1, 2)  + VV(I,2) 

for  each  unit  I = 1,2,...,  NU,  where  NU  is  the  total  number  of  units. 
See  Appendix  A for  a detailed  explanation  of  TP  and  W. 

After  each  unit  is  moved,  a check  is  made  to  determine  whether  it 
is  time  for  that  unit  to  change  its  direction  of  motion  (change  times  are 
stored  in  W(I,.3)).  If  so,  the  new  velocity  vector  is  retrieved  from  the 
matrix  W1  where  all  velocity  vectors  are  stored,  and  replaces  the  existing 
one  in  VV. 

The  second  stage  of  the  movement  computation,  that  of  computing  z, 
has  been  implemented  in  two  different  ways  in  WHERSM,  although  only  one 
appears  in  appendix  B.  (The  other  is  too  specific  to  include  here  and, 
furthermore,  is  trivial  to  implement.)  The  first  one,  the  one  not 
included  here,  is  simply  a matter  of  coding  up  a continuous  function 


) 

representation  for  the  z coordinate  and  computing  a value  of  z 
once  x and  y have  been  determined. 

The  method  of  determining  z that  is_  included  in  Appendix  B is 

somewhat  complicated,  in  that  it  involves  a preprocessor,  overlays,  a 

mass  storage  retrieval  routine,  and  parallel  processing.  The  preprocessor 

is  a routine  called  BUILD,  whose  function  is  to  read  a tape  containing 

topographical  data  in  the  format  designed  by  the  Electromagnetic  Compati- 

((*) **) 

bility  Analysis  Center  , decode  and  unpack  the  individual  altitudes 
above  sea  level  that  are  given  on  that  tape,  write  the  individual 
elevations  out  on  the  drum,  create  an  index  (NDX)  to  the  information  on 
the  drum,  and  initiate  the  overlaying  of  itself  by  the  rest  of  the  WHERSM 
code . 

The  subroutine  that  actually  provides  the  z coordinate,  HEIGHT, 
is  mainly  a drum  retrieval  routine.  HEIGHT  takes  the  x and  y that 
are  input  to  it,  computes  the  nearest  available  grid  point  to  that  x 
and  y,  determines  the  position  on  the  drum  of  the  altitude  associated 
with  that  grid  point,  and  initiates  the  retrieval  of  that  altitude.  This 
last  point  is  important  — HEIGHT  only  initiates  the  read.  It  does  not 
wait  for  the  completion  of  that  operation,  but  instead  immediately  returns 
control  to  the  calling  routine,  which  has  the  option  of  waiting  for  the 
completion  of  the  retrieval  or  of  performing  other  tasks  that  are  not 
dependent  on  knowledge  of  the  new  z coordinate.  In  view  of  the  fact 
that  the  time  required  to  retrieve  one  elevation  from  the  drum  is  approxi- 
mately 4.254  milliseconds,  the  above  parallel  processing  capability 
provided  by  the  HEIGHT  routine  allows  considerable  time  savings  to  be 


(*)  Note  that  obtaining  such  a function  for  "f ictitious"  terrains  is 
rather  easy,  but  that  surface  fitting  for  "real"  terrains  may  not 
be  trivial.  See,  e.g.,  Pavilidis,  T.,  "Piecewise  Approximations 
of  Functions  of  Two  Variables  and  its  Application  in  Topographic 
Data  Reduction",  Princeton  University,  Tech.  Report  No.  86,  Sept. 
1970,  and  "A  Program  Implementing  Piecewise  Linear  Approximations 
of  Functions  of  Two  Variables",  Princeton  University  Tech.  Memo. 
No.  3,  March  1971. 

(**)  Encoded  altitudes  above  sea  level  at  points  that  are  three  seconds 
of  latitude  and  longitude  apart. 
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achieved  through  the  thoughtful  organization  of  the  main  routine  relative 
to  height  retrieval. 

The  previous  comment  regarding  organization  serves  to  introduce  the 
last  point  to  be  discussed  in  this  section  — the  reason  for  including 
two  movement  modules  — one  (Box  2)  for  the  current  locatee,  and-  one 
(Box  3)  for  all  other  units  in  the  field.  Although  all  units  move  every 
subcycle,  a new  height  is  determined  only  when  a unit  serves  as  locatee. 

(In  between,  heights  are  treated  as  remaining  constant.)  Hence,  while 
the  height  for  the  locatee  is  being  retrieved  from  the  drum,  the  program 
can  continue  on  with  the  moving  of  all  other  units,  thus  taking  advantage 
of  the  parallel  processing  provided  by  HEIGHT.  There  is  very  little  loss 
of  "realism"  with  this  approach,  since  due  to  the  design  criteria  no  unit 
moves  farther  than  44  meters  between  times  when  it  is  a locatee,  and  further, 
since  the  topographical  grid  points  define  a rectangle  with  dimensions 
approximately  68m  x 92m.  Hence,  only  a small  distortion  of  an  already 
grossly  simplified  terrain  structure  occurs. 

There  is  yet  another  reason  for  separating  movement  into  two  modules. 
Namely,  it  reduces  the  effort  involved  in  implementing  a very  specific 
type  of  movement  scenario.  For  example,  the  code  in  Appendix  B includes 
a technique  for  keeping  units  from  spreading  out  too  much,  by  disallowing 
movement  past  a boundary  that  moves  at  a speed  determined  by  input.  (This 
is  convenient  when  a randomly  generated  scenario  that  includes  foot  soldiers 
and  ground  vehicles  is  used.  By  letting  the  boundary  move  at  foot  soldier 
speed,  the  ground  vehicles  can  be  kept  "close  by".)  In  this  case,  it  is 
wasteful  to  monitor  a unit’s  movements  all  the  time.  To  check  only  when 
a unit  serves  as  locatee  is  sufficient.  And  this  of  course  is  easier  and 
faster  when  there  is  a separate  movement  module  for  the  locatee. 


3.  BOX  4:  THE  RANGE  COMPUTATION 


The  range  computation  module  is  responsible  for  computing  reported 
ranges  by  computing  true  ranges  and  perturbing  them,  and  for  determining 
whether  each  range  should  be  reported.  The  computation  of  true  ranges 
involves  only  Euclidean  distance  calculations,  and  the  perturbation  values 
are  obtained  by  a call  to  subroutine  RG.  Subroutine  RG  generates  pseudo- 
random numbers  having  a Gaussian  distribution  with  mean  zero  and  variance 
2 

36m  , and  then  randomly  chooses  1%  of  them  to  which  it  imputes  a.  positive 

bias  by  adding  7 meters  to  each.  Except  for  the  biasing  procedure,  the 

method  used  to  generate  these  numbers  is  the  same  as  the  one  proposed  by 

(*) 

A.  Rotenburg  . RG  calls  RGU  which  provides  uniformly  distributed  pseudo 
random  numbers,  generated  by  the  multiplicative  congruential  method.  Qther 
distributions,  or  other  parameter  values  in  these,  could  easily  be  sub- 
stituted in  these  routines.  In  order  to  save  time, since  RGU  is  called  so 
frequently  during  the  simulation,  Rit  is  written  in  SLEUTH,  the  machine 
language  for  the  UNIVAC  1108.  However,  it  is  very  short  and  can  easily 
be  rewritten  in  FORTRAN. 

There  are  two  items  that  determine  whether  a particular  range  should 
be  reportad.  Tha  first  of  these  is  the  logical  matrix  RPT,  which  is 
initialized  at  input  time,  and  which,  for  each  unit,  has  the  logical 
value  ’’TRUE"  if  that  unit  is  one  that  reports  received  ranges,  and  "FALSE" 
if  it  is  not.  Note  that  the  maintenance  of  such  a matrix  allows  the  test- 
ing of  scenarios  in  which  the  set  of  reporting  units  changes  during  a run, 
and  where  such  changes  are  either  exogenous  jor  endogenous  events . 

The  second  item  that  determines  whether  a particular  range  should  he 
reported  is  the  matrix  of  intervisibilities  (IT) , which  records  whether 
or  not  a radio  line  of  sight  (LOS)  exists  between  each  pair  of  units 


(*)  Rotenburg,  A.,  "A  New  Pseudo-Random  Number  GeneratQr",  J.  ACM  7_. 
(1960),  pp.  75-77. 


identified  by  the  subscripts  for  IT.  Also,  each,  row  of  IT  is  condensed 
so  that  36  row  entries  fit  into  one  computer  word,  thereby  greatly  reducing 
storage  use.  Therefore,  the  only  concern  of  the  range  computation  module 
regarding  intervisibilities  is  that  of  locating  and  retrieving  a particular 
entry  in  IT  from  a block  of  36  entries,  which  informs  it  whether  the  range 
between  the  defining  units  should  be  calculated. 

In  summary,  then,  for  each  subcycle  the  range  computation  module 
steps  through  a list  of  the  units  (excluding,  of  course,  the  one  cur- 
rently being  located)  and  checks  RPT  to  see  first  if  a particular  unit 
is  a reporter.  If  so,  IT  is  checked  to  determine  whether  that  unit  can 
"see"  the  locatee.  If  that  also  is  true,  then  the  true  distance  is 
calculated,  RG  is  called  to  obtain  a perturbation  value,  and  that  value 
is  added  to  the  true  distance  yielding  a reported  range.  Reported  ranges 
are  stored  in  the  matrix  RANGE,  with  the  unit  numbers  of  the  units  that 
reported  each  range  in  matrix  LOCATR. 
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4.  BOXES  5,  6,  and  7;  THE  LOCATION  ATTEMPT,  STATISTICS 
ACCUMULATION,  AND  "END  OF  CYCLE"  TEST 


The  discussions  of  Boxes  5,  6,  and  7 have  been  combined  into  this 
one  section  simply  because  there  is  very  little  to  say  about  any  of  them. 

The  location  attempt  consists  of  a call  to  subroutine  FINDIT,  which 
is  a preprocessor  for  the  actual  location  algorithms.  As  was  mentioned 
in  Section  1,  no  discussion  of  the  algorithms  is  provided  here;  never- 
theless, some  discussion  of  the  input  to  the  algorithms  from  WHERSM 
should  be  given.  All  input  to  the  position  location  algorithms  is 
achieved  via  the  labeled  common  block  CB1,  which  consists  of  the  follow- 
ing variables:  EP,  LOCATR,  RANGE,  S,  TLOC,  W,  IU,  CLK,  XNEW,  YNEW,  ZNEW, 

FXD,  RPT,  and  NTYPE.  For  a description  of  each  of  these  variables,  see 
Appendix  A. 

The  error  statistics  that  are  accumulated  are:  the  current  coor- 

dinate errors  (coordinate-wise  subtraction  of  true  from  estimated  posi- 
tions) , the  locatee’s  position  location  error  (Euclidean  distance  between 
true  and  estimated  positions  in  xy-plane) , the  error  in  the  z coor- 
dinate, the  max  error  in  the  xy-plane  over  all  cycles  and  over  the  current 
block  of  KIC  cycles  (KIC  is  an  input  parameter) , the  maximum  error  in 
z over  all  cycles  and  over  the  current  KIC  cycles , and  the  time  to  per- 
form the  location  operation  (obtained  by  interrogating  the  computer’s 

internal  clock)  . The  statistics  accumulated  hare  are  used  to  produce 

the  intermediate  output  explained  in  Section  11. 

The  "end  of  cycle"  test  CBox  7 of  the  flowchart  In  figure  1)  is 
nothing  more  than  the  end  of  a DO  loop  which  steps  the  simulation  program 
through  each  one  of  the  subcycles  in  a cycle. 
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5.  BOXES  8 and  9:  THE" INTERVISIBILITY  TABLE  UPDATE" 


The  "intervisibility  update"  module  is  concerned  with  maintaining  an 
LOS  pattern  that  is  checked  by  the  WHERSM  program  during  range  computation 
(see  Section  3).  The  "current"  intervisibility  pattern  is  kept  in  the 
matrix  IT,  which  is  automatically  updated  in  a manner  determined  by  param- 
eters that  are  input  to  the  simulation  program.  This  section  explains  the 
manner  in  which  IT  is  updated. 

The  intervisibility  snapshot  matrices  that  have  been  discussed  in 
Sections  1 and  3 are  read  in  from  tape  at  input  time  and  stored  in  the 
matrix  ITB.  (The  format  of  the  tape  is  explained  in  Section  10.)  At 
the  start  of  a simulation  run,  IT  is  initialized  to  the  first  of  the 
snapshot  matrices — the  matrix  that  corresponds  to  a line  of  sight  snap- 
shot taken  at  time  zero.  Every  K2C  cycles  thereafter,  K3C  randomly 
chosen  rows  of  IT  are  overlaid  with  the  corresponding  rows  of  the  snap- 
shot matrix  that  is  next  in  sequence.  (K2C  and  K3C  are  input  parameters; 
see  Section  10 ♦)  This  overlaying  process  continues  until  simulated  time 
has  been  advanced  to  the  time  at  which  that  next  snapshot  was  taken,  at 
which  point  the  intervisibility  update  module  completes  the  overlaying 
of  all  rows  not  previously  updated,  and  prepares  to  start  over  again 
with  the  new  matrix  and  the  one  after  that. 

The  choice  of  a row  to  be  updated  is  achieved  by  choosing  a pseudo- 
random number,  call  it  I,  between  1 and  NU,  where  NU  is  the  number  of 
units  in  the  field  and  also  the  number  of  rows  and  columns  in  IT.  If  the 
corresponding  row  of  IT  has  not  previously  been  updated  (determined  by 
a check  of  the  logical  matrix  ITIND) , it  is  then  updated.  If  that  row 
has  been  updated,  try  row  1+1.  If  that  row  has  been  updated,  try  row 
1+2.  This  search  process  continues  until  either  an  "updatable"  row  is 
found  or  until  the  checking  process  circles  back  to  row  I.  (This  last 
comment  implies,  of  course,  that  row  NU+1  is  considered  to  be  row  1.) 

If  no  updatable  row  is  found,  nothing  changes  and  the  "intervisibility 
update"  module  has  completed  its  task. 
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6.  BOXES  10  and  11:  THE  PRODUCTION  OF  INTERMEDIATE  OUTPUT 


This  section  of  code  controls  the  production  of  intermediate  error 
output  as  portrayed  in  figure  3,  and  which  is  discussed  in  Section  11. 

Basically,  this  module’s  function  is  to  convert  the  statistics 
accumulated  in  Box  6 to  a form  that  is  more  meaningful  for  output.  For 
example,  sums  of  previous  errors  are  divided  by  the  appropriate  number 
to  obtain  averages,  coordinate  error  counts  are  made  to  produce  quadrant 
counts,  etc. 

The  frequency  with  which  the  intermediate  error  output  is  produced 
is  controlled  by  an  input  parameter  called  K1C.  The  output  summary  is 
provided  at  the  end  of  every  K1C  cycles. 


7.  BOX  12:  UPDATE  THE  STATIONARY  UNIT  LIST 


The  stationary  unit  list  (logical  matrix  FXD)  is  a list  that  contains 
''known"  information  regarding  which  units  are  moving  and  which  units  are 
not  moving  throughout  the  simulation.  By  "known"  is  meant  that  this 
information  i_s_  available  to  the  position  location  algorithms. 

There  are  two  types  of  real  world  scenarios  that  the  matrix  FXD  has 
been  instrumental  in  simulating.  The  first  of  these  is  the  situation 
where  units  stop  and  start  throughout  the  "battle",  and  send  a signal 
back  to  the  master  computer  indicating  that  they  have  just  stopped  or 
that  they  have  just  started  moving.  The  second  scenario  is  one  in  which 
certain  units  are  "surveyed  in"  and  never  move,  whereas  no  other  units 
are  ever,  treated  as  fixed.  Knowledge  of  these  facts  is  very  helpful  to 
the  position  location  algorithms , but  the  "how"  and  the  "why"  of  that 
helpfulness  will  not  be  discussed  here.  For  that  information,  the  reader 

is  referred  to  Section  3 of  the  parent  document  to  which,  this  Working 
Paper  is  an  appendix. 

In  the  current  version  of  WHERSM,  the  starting  and  stopping  of  units 
mentioned  above  is  achieved  by  maintaining  the  matrix  FXD,  and  by  updating 
it  in  a random  manner  that  is  controlled  by  parameters  which  are  input  to 
the  simulation.  For  each  unit,  FXD  is  "TRUE"  if  that  unit  is  currently 
fixed,  and  "FALSE’'  if  that  unit  is  currently  moving.  The  input  parameter 
PSS(2)  is  the  average  time  a unit  remains  fixed  once  it  stops  moving, 
while  PSS(l)  is  the  average  time  it  keeps  moving  once  it  starts  moving. 
These  times  are  converted  to  probabilities,  stored  back  into  PSS(l)  and 
PSS(2),  and  those  probabilities  in  turn  are  used  to  determine,  at  each 
cycle,  whether  a unit  should  have  its  state  ("fixed"  or  "moving")  changed. 

The  determination  of  state  changes  is  precisely  the  purpose  of  the 
"stationary  list  update"  module.  It  is  achieved  by  obtaining  a uniformly 
distributed  (0,1)  pseudo-random  number  from  the  internal  subroutine  RU, 
and,  if  that  number  is  less  than  the  appropriate  PSS  probability  value, 
changing  the  state  of  that  particular  unit. 
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It  should  be  noted  that  in  the  second  scenario  mentioned  above, 
one  of  course  does  not  want  the  stationary  unit  list  to  be  changed  at 
all  during  the  simulation.  This  can  be  guaranteed  simply  by  letting 
the  values  of  PSS(l)  and  PSS(2)  at  input  time  be  very  large  numbers. 
It  should  also  be  noted  that  FXD  is.  initialized  at  input  time. 


8.  BOXES  13  and  14:  PRODUCE  THE  X-Y  PLOTS 

This  module  is  very  small,  consisting  mainly  of  two  calls  to  sub- 
routine GRAPH.  The  first  call  produces'  x-y  plots  of  the  true  positions 
of  all  the  units,  and  the  second  call  produces  x-y  plots  of  the  esti- 
mated positions  of  all  the  units. 

There  are  eight  parameters  in  the  call  statement  for  subroutine 
GRAPH:  N,  X,  Y,  IWRD,  XMAX,  XMIN,  YMAX,  and  YMIN.  N is  the  number  of 
points  to  be  plotted,  X and  Y are  vectors  of  the  x and  y coor- 
dinates of  the  points  to  be  plotted,  IWRD  contains  either  ’TRUE’  or  ’EST.' 
depending  on  whether  the  points  to  be  plotted  are  true  or  estimated  posi- 
tions, and  XMAX,  XMIN,  YMAX,  and  YMIN  are  the  upper  and  lower  limits  of 
the  axes.  Those  last  (XMAX  through  YMIN)  are  part  of  the  input  to  the 
simulation  program. 

For  a description  of  the  plots  that  are  produced,  see  Section  11. 

The  frequency  with  which  the  plots  are  produced  is  controlled  by 
the  input  parameter  K4C ; the  plots  are  produced  after  every  K4C  cycles. 
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9.  BOX  15;  THE  END  OF  SIMULATION 


Box  15  of  the  flowchart  in  figure  1 is  the  only  one  that  does  not 
correspond  to  an  easily  recognizable  module  in  the  WHERSM  code.  Rather, 
its  function  is  performed  by  a set  of  IF  statements  scattered  throughout 
the  other  modules.  The  IF  statements  are  all  identical;  they  test 
whether  the  number  of  cycles  required  to  run  the  simulation  for  the 
length  of  time  indicated  by  the  input  parameter  SIMTIM  have  been  completed. 
If  so,  then  a final  set  of  "intermediate"  output  is  produced  before 
termination. 


10.  THE  INPUT  TO  THE  SIMULATION 


There  are  three 

types  of  input  used 

in  the 

WHERSM 

program:  input 

from  cards,  tape  input,  and  "hardwired" 

input  in  the  form  of  DATA  and 

PARAMETER  statements 

inserted  in  the  program  deck.  This  section  will 

discuss  each  type  of 

input  with  regard  to  input 

values 

and  formats. 

Card  input  is  explained  first. 

The  "Catch-All" 

Cards : These  cards 

contain  the  various  parameters 

and  constants  needed  by  the  program*  there  are 

two  of 

these  cards , and 

the  table  below  explains  their  contents. 

For  a 

definition  of  the  vari- 

ables  mentioned,  see 

Appendix  A. 

Columns 

Format 

Variable  Name 

Card  No.  1: 

1-5 

15 

NU 

6-10 

15 

NS 

11-15 

15 

K1C 

16-20 

15 

K2C 

21-25 

15 

K3C 

26-30 

15 

K4C 

31-35 

15 

INR  . 

36-40 

. 15 

NSP 

41-50 

F10.5 

PSS(l) 

51-60 

F10.5 

PSS  (2) 

61-70 

F10.5 

BND 

71-80 

F10.5 

BNDSP 

Card  No.  2: 

. 1-10 

F10.5 

CYCTM 

11-20 

F10.5 

SIMTIM 

21-30 

F10.5 

> 

XMAX 

31-40 

F10.5 

XMIN 

41-50 

F10.5 

YMAX 

51-60 

F10.5 

YMIN 

The  "Reporting  Units"  Cards:  These  cards  contain  the  unit  numbers 

of  the  units  that  will  be  used  at  the  start  of  the  simulation  as  reporters 
More  specifically,  the  elements  of  the  matrix  RPT  that  correspond  to  the 
numbers  that  appear  on  the  Reporting  Units  Cards  are  the  only  elements 
that  are  initialized  to  the  logical  value  "TRUE".  All  others,  are'  ini- 
tialized to  "FALSE".  There  are  20  numbers  per  card,  each  with  a format 
of  14.  The  number  of  unit-numbers  that  appear  on  these  cards  is  the  value 
of  INR  (see  Catch-all  cards  above).  Consequently,  there  will  be 
[(INR-l)/20]  + 1 Reporting  Units  cards. 

The  "Fixed  Units"  Cards:  These  cards  contain  the  unit  numbers  of 

the  units  that  will  be  used  at  the  start  of  the  simulation  as  fixed  units. 
The  elements  of  the  matrix  FXD  are  initialized  in  the  same  manner  as  are 
the  elements  of  the  matrix  RPT  above.  The  format  for  each  card  again  is 
2014,  and  there  will  be  [(NSP~l)/20]  +1  of  these  cards. 

The  "Type"  Cards:  The  type  cards  contain  the  type  number  for  each 

unit  in  the  simulation.  There  currently  are  four  unit  types: 


Type  Number 


TyPe 


Speed 


1 

2 

3 

4 


Foot  Soldier 
Ground  Vehicle 
Low  Performance  Aircraft 
High  Performance  Aircraft 


11/2-2  km/h 
8-10  km/h 
161  km/h 
322  km/h 


There  will  be  exactly  NU  of  these  type  numbers  each  with  a format  of  14, 
and  therefore  [ (NU  - l)/20]  + 1 Type  cards. 

The  "Movement"  Cards:  These  cards  contain  all  the  information  needed 

by  the  program  so  that  it  can  provide  for  the  movement  of  units.  There 
will  be  one  of  these  cards  for  each  unit  in  the  scenario,  hence  NU  cards 
in  all.  Their  contents  are  as  follows: 


Card  Columns 

Format 

Meaning 

1-3 

13 

Unit  number 

4-21 

3A6 

Unit  ID  (anything) 

Starting  latitude: 

22-24 

13 

Degrees 

25-26 

12 

Minutes 

27-28 

12 

Seconds 

Starting  longitude: 

29-31 

13 

Degrees 

32-33 

12 

Minutes 

34-35 

12 

Seconds 

37-40 

F4.0 

Time  ^ ^ (hrs . ) 

41-44 

F4.0 

Azimuth  (deg.) 

45-48 

F4.0 

Speed  (km/h) 

50-53 

F4.0  • 

(it) 

Time  (hrs.) 

54-57 

F4.0 

Azimuth  (deg.) 

58-61 

F4.0. 

Speed  (km/h) 

63-66 

F4.0 

(it) 

Time  (hrs . ) 

67-70 

F4.0 

Azimuth  (deg . ) 

71-74 

F4.0 

Speed  (km/h) 

There  are  two  kinds  of  tape  input  used  in  the  WHERSM  program.  The 
first  of  these  is  the  tape  that  contains  the  height  information  in  the 
form  of  a digitized  terrain  map.  The  format  of  this  tape  is  as  appears 
in  figure  2a  but  no  further  discussion  of  format  will  be  given  here  since 
the.  terrain  tape  was  obtained  directly  from  the  Electromagnetic  Compatability (*) 


(*)  By  "time"  is  meant  the  length  of  time  during  which  the  immediately 
following  azimuth  and  speed  are  used.  Control  is  passed  to  the 
next  set  when  "time  runs  out"  for  this  set. 
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(*) 

Analysis  Center  which  has  published  a report  that  documents  such  tapes 
in  detail.  The  terrain  tape  must  be  mounted  on  logical  unit  number  8. 

The  other  tape  used  as  input  to  WHERSM  is  one  that  is  mounted  on 
logical  unit  number  7 and  contains  the  intervisibility  snapshot  matrices 
discussed  in  Sections  1 and  3.  It  is  a binary  tape,  created  as  follows: 
Let  ITB  be  a three-dimensional  matrix  that  contains  the  set  of  "packed” 
intervisibility  snapshot  matrices; 

ITB (I , J ,1) , I = 1,2,. ..,NU;  J = 1,2,...,  [ (NU  - l)/36]  + 1 

is  the  snapshot  matrix  taken  at  the  first  discrete,  time  value;  ITB(I,J,2) 
is  the  matrix  taken  at  the  second  discrete  time  value;  etc.  Then  the 
tape  should  have  been  created  using  a write  statement  of  the  form; 


WRITE  (7)  ( ( (ITB (I , J ,K) , J = 1,  NITCOL) , I = 1,NU),  K = 1,K1) 

where  NITCOL  = [ (NU  - 1 ) / 36 ] + 1 (implying  of  course  that  36  entries 
of  each  row  of  the  2 dimensional  matrices  are  packed  into  one  computer 
word),  and  K1  equals  the  number  of  snapshot  matrices  that  exist. 

The  third  and  last  type  of  input  required  to  make  a simulation  run 
is  the  type  that  must  physically  be  inserted  in  the  program  deck  - the 
PARAMETER  and  DATA  statement  cards . 

The  table  below  lists  these.  The  reader  is  referred  to  Appendix  A 
for  their  meanings. 


(*)  Cleramitt,  M.  and  Stone,  R. , "The  Topographic  Information  System", 

Tech.  Note  No.  ECAC-TN-007-222 , June  1969,  Electromagnetic  Compati- 
bility Analysis  Center,  Annapolis,  Md. 


f ' , 


Type  of  Statement  Routine 


Variable  Name 

Contained  In 

Contained  In 

MNLAT 

DATA 

BUILD 

MXLON 

DATA 

BUILD 

TPOUT 

DATA 

WHERSM 

IFREQ 

DATA 

WHERSM 

CNX 

DATA 

HEIGHT 

CNY 

DATA 

HEIGHT 

W 

DATA 

WHERSM 

FXDINI 

DATA 

WHERSM 

LI 

PARAMETER 

WHERSM 

L2 

PARAMETER 

WHERSM 

L6 

PARAMETER 

BUILD 

WHERSM 

HEIGHT 

L7 

PARAMETER 

BUILD 

WHERSM 

HEIGHT 

11.  THE  OUTPUT  FROM  THE  SIMULATION 


Figure  3 contains  a sample  of  the  intermediate  output  produced  by 
WHEPvSM.  This  section  explains  the  contents  of  that  output  sample. 

The  first  line  of  output  contains  the  cycle  number  of  that  cycle 
after  which  the  summary  was  produced,  the  simulation  clock  time,  the 
average  real  clock  time  the  program  took  to  provide  for  the  processing 
of  one  sub cycle,  and  the  average  time  taken  by  the  algorithms  to  perform 
one  location  operation. 

Each  row  of  the  table  contains  statistics  for  the  unit  whose  unit 
number  appears  in  column  1.  Columns  2-7  contain  statistics  accumulated 
over  all  cycles  previous  to  the  current  one.  They  are, in  order: 

the  average  error  in  the  xy-plane, 
the  average  error  in  the  z coordinate, 
the  maximum  error  in  the  xy-plane, 

the  cycle  during  which  that  maximum  error  occurred, 

the  maximum  z error  that  occurred,  and 

the  cycle  during  which  that  maximum  z error  occurred. 

Columns  8-11  contain  statistics  accumulated  over  the  previous  K1C 
cycles.  Those  statistics  are:  the  average  x-y  error,  the  average  z 

error,  the  maximum  x-y  error  that  occurred,  and  the  cycle  during  which 
that  error  occurred. 

Columns  12-14  contain  the  individual  coordinate  errors  (estimated 
minus  true)  that  occurred  the  last  time  each  unit  was  located  prior  to 
the  printing  of  the  intermediate  output.  Columns  15,  16,  and  17  list 
the  current  values  of  the  S,RPT,  and  FXD  vectors  respectively;  and 
column  18  gives  the  current  value  of  the  W vector,  the  weight. 

The  last  row  of  the  table,  which  begins  with  the  word  "ALL",  gives 
for  a column  of  averages  the  average  of  the  numbers  appearing  above  it; 
and  for  a column  of  maximums,  the  maximum  of  the  maximums  in  the  column 
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The  next  item  of  output  is  the  quadrant  counts  of  errors.  This  is 
a display  of  the  directions  of  errors.  It  is  obtained  from  columns  12  and 
13  of  the  table  by  counting  the  number  of  entries  in  those  columns  (by 
their  algebraic  signs)  that  fall  in’  each  quadrant. 

The  last  items  of  output  are  the  xy-plots,  an  example  of  which  is 
given  in  figure  4.  The  two  lines  at  the  top  of  the  graph  state  whether 
the  graph  is  of  true  or  estimated  positions,  list  the  unit  numbers  of 
the  units  that  appear  on  that  graph,  and  list  the  plotting  symbols  used 
to  represent  those  units.  The  letters  of  the  alphabet  are  used  as  plotting 
symbols,  so  that  if  more  than  26  units  are  in  the  simulation,  a separate 
graph  will  be  produced  for  every  26  (or  fraction  thereof)  units.  If 
more  than  one  unit  occupies  the  same  graph  position,  a number  will  appear 
in  that  graph  position  indicating  the  number  of  units  in  that  position. 

The  axis  ranges  are  controlled  by  XMAX,  }£MIN,  YMAX,  and  YMIN,  which  are 
part  of  the  input  to  the  simulation,  and  if  points  fall  beyond,  those 
boundaries  an  annotation  to  the  graph  will  appear  indicating  the  occur- 
ence of  that  situation. 
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Figure  4:  Example  of  XY-Plots 


APPENDIX  A:  A DICTIONARY  OF  VARIABLES  USED  IN  THE  SIMULATION  PROGRAMS 


In  this  appendix,  definitions  of  all  non-trivial  variables  used  in 
the  three  main  simulation  routines  (BUILD,  WHERSM,  HEIGHT)  are  given. 

In  each  case,. the  variable  being  defined  is  underlined,  and  if  it  is  a 
dimensioned  variable,  it  is  written  exactly  as  that  variable  should  appear 
in  a DIMENSION  statement.  Immediately  following  the  variable  name  will 
appear  either  a B,W,  or  H (or  any  combination  thereof).  By  this  is 
meant  that  the  variable  being  defined  is  a variable  of  BUILD,  WHERSM 
or  HEIGHT. 

Note:  All  variables  conform  to  the  FORTRAN  implicit  naming  con- 

ventions regarding  mode. 

***** 

AVASBT  W.  The  sum  of  real  times  used  to  perform  position  estimations. 

AVKTR  W.  The  number  of  time  intervals  accumulated  in  AVASBT  and  AVSUBT. 
AVSUBT  W.  The  sum  of  real  times  used  to  simulate  a "real  world"  sub- 
cycle (includes  AVASBT) . 

AZ  (4)  W.  Buffer  used  in  reading  in  movement  azimuths. 

BNP  W.  Boundary  beyond  which  units  are  not  allowed  to  move.  Unit  of  mea- 
surement is  the  kilometer.  If  this  option  is  not  desired,  set  BND  to 
some  large  number. 

BNPSP  W.  On  input,  the  speed  (km/h)  with  which  BND  moves  in  both  directions 
Immediately  converted  to  a value  that  is  added  to  BND  at  every  sub cycle. 

0(14641)  B.  An  output  buffer  that  contains  decoded  elevations  above 


sea  level. 


■ 


CLK  W.  The  simulation’s  clock.  Contains  simulated  time  in  seconds. 

CNX  H.  The  inverse  of  the  average  longitudinal  distance  (in  meters) 
between  grid  points  on  the  terrain  tape. 

CNY  H.  Same  as  CNX  but  for  latitude, 

CT (4)  W.  Input  buffer  used  to  read  in  the  length  of  time  a unit  is 
to  move  according  to  a particular  azimuth  and  speed. 

CYCTM  W.  The  "real  time"  length  of  a cycle  in  seconds. 

D (14641)  B.  Output  buffer  that  contains  decoded  altitudes  above-  sea 
level . 

DELTA  W.  "Real  time"  length  of  a subcycle. 

DX(L1)  W.  List  used  to  store  the  difference  in  the  x-coardlnates  of 
the  true  and  estimated  positions.  Computed  immediately  after  a unit 
has  had  its  position  estimated. 

DY (LI)  W.  Same  as  DX,  but  for  y-coordinates . 

DZ  (LI)  W.  Same  as  DX,  but  for  z-coordinates . 

EP (LI ,14)  W.  Contains  estimated  positions  for  each  unit,  in  addition 
to  some  data  relevant  to  those  positions.  EP(I,l-3)  is  the  estimated 
position  of  unit  I which  was  determined  the  next-to-last  time  it  served 
as  locatee;  EP(I,4-6),  the  last  time  it  served  as  locatee.  The  rest 
of  the  EP  matrix  is  used  to  manipulate  positions  when  the  relocation 
option  is  used  (see  Working  Paper  No.  8). 

ERR  W.  The  error  in  position  estimation  in  the  xy-plane  and  then  in 
the  z-coordinate.  The  Euclidean  distance  between  true  (TP)  and  estimated 
(EP)  positions  is  used  as  the  error. 


" 


IFY  H.  Same  as  IFX,  but  for  the  y-coordinate.  Units  of  measurement 
here  are  3 seconds  of  latitutde. 

IGP  W.  Pointer  to  the  unit  with  the  largest  xy-error  over  all  cycles. 

IGT  W.  Pionter  to  the  unit  with  the  largest  xy-error  over  the  last  K1C 
cycles . 

INR  W.  Initial  number  of  units  allowed  to  report. 

IP  W.  Pointer  used  in  constructing  Wl. 

IPX  H.  The  x-coordinate  of  some  other  position  (for  which  a height  is 
already  available) , after  it  has  been  converted  to  the  units  of  the 
terrain. 

IPY  H.  Same  as  IPX,  but  for  y-coordinate. 

IQO  W,  Number  of  units  whose  estimated  position  is  exact. 

IQ1  W.  The  number  of  units  whose  estimated  position  is  in  error  In 
the  direction  of  the  first  quadrant. 

IQ2  W.  Same  as  IQ1,  but  for  quadrant  2. 

IQ3  W.  Same  as  IQ1,  but  for  quadrant  3. 

IQ4  W.  same  as  IQ1,  but  for  quadrant  4. 

IQUAN  B.  Quantization  factor  for  the  block  of  terrain  grid  points 
currently  being  decoded.  (See  ECAC^TNt'007t-222 , op . cit . ) 

IR  W.  Counter  for  the  number  of  ranges  computed  in  a particular  sub cycle 
IT (L1,L3)  W.  The  intervisibility  table  currently- being  used. 

ITB(L1 ,L3 ,7)  W.  The  set  of  all  intervisibility  tables. 

ITIND (LI)  W.  Logical  matrix  which  indicates  for  each  row  of  IT,  whether 


that  row  has  been  updated  yet. 


IU  W.  The  unit  number  of  the  unit  currently  serving  as  locatee. 

IV  W.  The  CO-1)  answer  to  the.  question  of  existence  of  an  LQS  between 
the  pair  of  units  in  question. 

IVMAT  W.  Pointer  to  the  in ter visibility'  matrix  currently  being  used  to 
update  IT . 

IZT  W.  Pointer  to  the  unit  with  the  largest  z-error  over  all  cycles. 
K1C  W.  The  number  of  cycles  after  which  intermediate  output  is  to  be 
produced. 

K1CC  W.  Counter  used  with  K1C. 

K2C  W.  The  number  of  cyeles  after  which  an  intervisibility  update 
is  to  occur. 

K2CC  W.  Counter  used  with  K2C. 

K3C  W.  The  number  of  rows  of  IT  to  be  updated  when  an  in ter visibility 
update  occurs. 

K4C  W.  The  number  of  cycles  after  which  xy  plots  are  to  be  produced. 
K4CC  W.  Counter  used  with  K4C. 

KUF (LI)  W.  For  I = 1,2,..., NS,  KUF(I)  is  the  unit  number  of  the 
unit  which  is  to  serve  as  locatee  during  the  Ith  subcycle. 

KZ  W.  The  seed  value  for  the  internal  uniform  psuedo -random  number 
generator,  RU. 

B.  Status  parameter  for  I/O  processor  NTRAN.  L equals^  -1  until 
the  completion  of  the  I/O  command. 

LI  W.  Parameter  variable  used  in  dimensioning.  LI  must  be  greater 
than  or  equal  to  NU. 


' 


■ 


L2  W.  Parameter  variable  us.ed  in  dimensioning.  L2  NS. 

L3-L5  W.  Parameter  variables  used  in  dimensioning,  but  since  they  are 
strictly  dependent  on  LI,  they  are  automatically  set. 

L6  B,W,H.  Conceptually  organize  all  the  6’  x 6’  terrain  blocks  on 
the  terrain  tape  according  to  the  latitude-longitude  of  their  southwest 
corner.  L6,  then,  is  the  number  of  these  blocks  when  counting  from 
west  to  east. 

L7  B,W,H.  Same  as  L6,  but  count  from  south,  to  north., 

LQCATR (LI)  W.  The  unit  numbers  of  the  units  whose  ranges  are  being 
reported  in  a given  subcycle. 

M W.  Status  parameter  for  I/O  processor.  See  L. 


Ml  B. 

Status 

parameter 

for 

I/O 

processor . 

See  L. 

M2  B. 

Status 

parameter 

for 

I/O 

processor . 

See  L. 

MAXCYC  W.  The  length  of  the  simulation  in  terms  of  cycles. 

MNLAT  B.  The  minimum  latitude  (in  seconds)  for  which,  a height  is 
available  on  the  terrain  tape. 

MXLON  B.  The  maximum  longitude  Cin  seconds)  for  which,  a height  is 
available  on  the  terrain  tape. 

NBITS  B.  .The  number  of  bits  into  which  the  elevation  in  question  is 
"packed". 

NCYC  W.  The  cycle  number  of  the  cycle  currently  being  processed. 


NDX(L6,L7)  B,H.  An  index  to  terrain  blocks  that  are  written  on  the 
drum  by  BUILD.  NDX(1,1)  is  the  position  of  the  6T  x 6*  terrain  block 
whose  southwest  corner  is  identified  by  (JINLAT,  MXLON) ; NDXClj2)‘  is  the 
position  of -block  (MNLAT,  MXLON  + 3600);  NDX(2,1)  is  block  (MNLAT  + 3600, 
MXLON);  etc. 

NITCOL  W.  The  number  of  computer  words  in  a packed  row  of  the  inter— 
visibility  matrices. 

NS  W.  The  number  of  subcycles  in  a cycle. 

NSP  W.  The  number  of  units  which,  at  the  start  of  the  simulation,  will 
be  identified  as  being  stationary. 

NTYPE (4)  W.  The  unit  type  of  each  unit.  If  NTYPE(I)  is  1,2,3,  or  4, 
then  unit  I is  respectively  a foot  soldier,  a ground  vehicle,  a low 
performance  aircraft,  or  a high  performance  aircraft. 

NU  W.  The  number  of  units  in  the  simulation. 

NUMEL  B.  The  number  of  elevations  per  word  for  a given  terrain  block 
represented  by  a logical  record  on  the  terrain  tape. 

PSS(2)  W.  At  input  time,  PSSCl).  is  the  average  time  a unit  keeps  moving 
once  it  starts  moving,  and  PSS(2)  is  the  average  time  it  remains  fixed 
once  it  stops.  Subsequently,  they  are  the  probabilities  of  changing 
state.  If  no  random  starting  and  stopping  is  desired,  set  botlrPSSCl) 
and  PSS(2)  to  very  large  numbers  on  the  catch-all-cards. 

RADS  W.  Constant  representing  the  number  of  degrees  in  a radian. 

RANGE (LI)  W. . The  ranges  being  reported  for  a given  subcycle. 


■ 


RCON  W.  Constant  used  in  converting  a speed  and  azimuth.  combination 


into  direction  cosines. 

RPT (LI)  W.  Logical  matrix  that  indicates  for  each  unit,  whether  or 
not  that  unit  is  a unit  that  is  allowed  to  report  ranges. 

S (LI)  W.  Logical  matrix  that  indicates  for  each  unit  whether  or  not 
that  unit  was  "successfully"  located  the  last  time  a position  estimation 
was  attempted.  S is  set  by  the  position  location  algorithms. 

SIMTIM  W.  The  "real  time"  length,  of  the  simulation,  in  minutes. 

SP (4)  W.  Input  buffer  used  in  reading  in  speeds  for  each  unit. 

STAT1 (LI ,10)  W.  For  each  row:  column  1 is  the  sum  of  position 

estimation  errors  (in  xy-plane)  over  all  cycles  for  the  unit  associated 
with  that  row;  column  2,  the  maximum  xy-error  that  occurred  over  all 
cycles  for  that  unit;  column  3,  the  cycle  number  of  the  cycle  during 
which  that  maximum  error  occurred;  4,  the  partial  sum  of  xy-errors 

over  the  current  KlC  cycles;  5,  the  maximum  xy-error  over  the  current 

* 

KlC  cycles;  6,  the  cycle  number  associated  with,  that  error;  7,  the 
sum  of  the  z-errors  over  all  cycles  for  that  unit;  8,  the  maximum 
z-error  over  all  cycles;  9,  the  cycle  number  associated  with  that 
z-error;  and  10,  the  partial  sum  of  z-errors  over  the  current  KlC 
cycles . 

STAT2 (LI, 2)  W.  Column  1 contains,  for  each  unit,  the  number  of  times 
that  unit  was  "successfully"  located  since  tire  start  of  the  simulation; 
and  column  2,  the  number  since  the  start  of  the  current  KlC  cycles. 


' 


TLOC (LI ,4)  W.  The  clock  time  at  which  the  estimated  positions,  in 
EP  were  obtained. 

TP  (LI ,3)  W.  The  current  true  position  of  each  unit.  TP  (1,1)  is  the 
x-coordinate;  TP  (1 ,2),  the  y-coordinate ; and  TP  (.1,3),  the  z-coordinate . 
TPOUT  W.  Logical  variable  that  allows  for  the  writing  of  all  inter- 
mediate output  on  logical  unit  9,  in  addition  to  the  standard  output 
unit. 

T1-T4  W.  Parameters  in  calls  to  the  system’s  clock.  Used  in  determining 
AVASBT  and  AVSUBT. 

W (LI ,4)  W.  Currently  used  velocity  vector  for  each  unit.  Column  1 
is  the  x component;  column  2,  the  y component;  3,  the  time  at 
which  that  vector  should  no  longer  be  used;  and  4,  a pointer  into  VV1 
where  the  next  velocity  vector  to  be  used  resides. 

VV1 (L5 ,4)  W.  The  set  of  all  velocity  vectors  for  all  units.  The 
format  of  each  vector  is  the  same  here  as  in  VV,  but  the  h (h  3 
currently)  vectors  for  unit  number  1 are  in  rows  1 through  h of  Wl, 
followed  by  the  vectors  for  unit  number  2,  etc. 

W(L1)  W.  Vector  used  by  the  position  location  algorithms  to  weight 
the  various  estimates  (see  Working  Papers  No.  7 and  8).  W should  be 
initialized  to  the  weight  to  be  assigned  to  non-fixed  units  at  the 
start  of  the  simulation. 

XMAX,  XHIN,  YMAX,  YMIN  W.  The  upper  and  lower  bounds  of  the  x-and 
y-axes  (in  meters)  to  be  used  in  producing  the  xy-^plots. 


. 

. 


. 


ZPSUM  W.  The  sum  of  the  position  estimation  errors  (in  the  z-coordinate) 


over  all  units  over  the  last  K1C  cycles.. 

ZTSUM  W.  The  sum  of  the  position  estimation  errors  (in  z-coordinate) 


over  all  units  over  all  cycles. 


v 


. 


APPENDIX  B:  LISTING  OF  THE  SIMULATION  PROGRAM 
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APPENDIX  C:  FLOWCHARTS  OF  THE  MAIN  SIMULATION  PROGRAMS 
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